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ABSTRACT 


This  volume  presents  a  study  conducted  to  determine  the  effects  thrust  bearings 
have  on  rotor-bearing  stability.  A  computer  program  was  written  in  order  to 
study  these  effects  and  permits  the  inclusion  of  the  thrust  bearing  character¬ 
istics  into  the  rotor  system.  A  manual  is  provided  for  the  program,  containing 
a  listing  of  the  program  and  detailed  instructions  for  preparation  of  input 
data.  A  technique  is  also  presented  which  permits  the  user  a  convenient  method 
of  evaluating  the  stability  of  a  rotor  system  with  and  without  the  thrust  bear¬ 
ing  data.  Extensive  design  data  are  presented  for  gas-lubricated,  externally 
pressurized  thrust  bearings. 
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SECTION  I 
INTRODUCTION 

For  rotors  possessing  a  high  degree  of  dissysasetry  or  for  rather  short  rotors  the 
Instability  whirl  motion  will  have  a  substantial  conical  whirl  component.  Under 
such  conditions  the  dynamic  forces  imposed  on  the  rotor  by  the  thrust  bearing  will 
have  an  Important  effect  on  the  dynamical  characteristics  and  the  stability  of  the 
rotor-bearing  system. 

The  purpose  of  the  study  described  in  this  report  was  to  develop  the  technology  fa 
determining  effects  of  thrust  bearing  forces  on  the  stability  of  rotors  supported 
in  fluid  film  bearings  and  to  assess  the  significance  of  this  effect  on  rotors  of 
practical  design.  To  this  end,  a  computer  program  was  developed  to  perform  many 
of  the  complex  calculations  necessary  for  determination  of  the  stability  limits  of 
complicated  rotor-bearing  systems  which  takes  account  of  effects  of  thrust  bearing 
stiffness  and  damping.  This  computer  program  is  an  extension  and  modification  of 
a  similar  program  developed  in  an  earlier  task  performed  under  contract  AF33(615) 
3238  (Ref.  1).  A  complete  description  of  the  program  including  a  computer  listing 
is  provided  in  the  present  report. 

The  present  report  provides  angular  stiffness  data  for  hydrostatic  thrust  bearing 
geometries  for  use  in  the  determination  of  rotor  stability.  Also  a  useful,  simpli 
fled  method  for  determining  stability  of  rotors  supported  on  bearings  all  of  which 
are  the  same  is  described.  Finally,  several  numerical  examples  of  the  calculation 
of  thrust  bearing  effects  on  rotor  stability  are  presented  and  the  significance  of 
these  effects  is  discussed. 
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SECTION  II 


THE  DETERMINATION  OF  THRESHOLD  SPEED  FOR  WHIRL  INSTABILITY  OF  A  ROTOR 

1.  General  Discussion 

As  is  well  established,  a  fluid  film  bearing  supporting  a  rotor  may  be  likened  to 
a  spring-dashpot  system  in  that  the  bearing  reaction  to  small  displacements  may 
be  expressed  in  terms  of  stiffness  and  damping  coefficients.  A  mass  supported  by 
springs  has  a  number  of  natural  spring-mass  frequencies  depending  on  the  complexity 
of  the  spring  system  and  the  possible  modes  of  motion  of  the  mass.  If  any  mode  of 
motion  is  excited  at  the  natural  frequency  by  an  external  harmonic  force,  then  the 
response  of  the  system  will  be  at  a  maximum  but,  if  the  system  possesses  effective 
damping  for  this  mode  of  excitation,  then  the  response  will  be  bounded  and  the  sys¬ 
tem  will  not  be  unstable.  If,  on  the  other  hand,  the  natural  frequency  of  a  mode 
of  motion  corresponds  to  a  condition  where  the  effective  damping  of  the  system  is 
negative,  then  the  motion  of  the  system  at  that  frequency  will  Increase  without 
bound  without  any  external  excitation,  and  the  system  is  considered  as  being  un¬ 
stable,  i.e.  subject  to  self-excited  vibration.  If  a  natural  frequency  occurs  at 
a  condition  of  zero  damping,  then  the  system  can  be  said  to  be  at  the  threshold  of 
instability,  where  Infinitesimally  small  exciting  forces  applied  over  a  period  of 
time  will  result  in  ever  increasing  amplitudes  of  response. 

The  so-called  critical  speeds  of  a  rotor-bearing  system  refer  to  resonant  responses 
of  the  system  to  unbalance  forces  of  the  system  which,  by  their  nature,  are  applied 
at  a  frequency  synchronous  with  the  rotational  speed  cf  the  shaft.  Since  fluid  film 
bearings  in  general  have  positive  damping  to  synchronous  speed  oscillations,  criti¬ 
cal  speeds  are  not  instabilities  but  simply  represent  a  condition  of  large  resonant 
response  to  inherent  unbalance  forces.  Fluid  film  bearings,  however,  do  have  an 
effective  damping  to  various  modes  of  motion  which  does  tend  toward  zero  when  the 
motion  occurs  at  some  fraction  of  the  running  speed,  usually  at  half  the  running 
speed.  Thus,  as  the  running  speed  of  a  rotor  increases  beyond  the  first  critical 
speed  or  natural  frequency  speed,  the  rotor  will  begin  to  approach  the  condition 
where  the  first  resonant  or  natural  frequency  of  the  system  will  coincide  with  the 
fractional  frequet  y  at  which  effective  damping  goes  to  zero.  When  this  coinci¬ 
dence  occurs,  the  rotor  is  said  to  have  reached  the  whirl  threshold  speed.  A 
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further  increase  in  rotor  speed  usually  results  in  a  very  rapid  growth  of  whirl 
amplitude  and,  in  almost  all  cases,  the  whirl  threshold  speed  represents  the 
upper  limit  for  safe  operation  of  the  rotor. 

He  might  note  at  this  time  that  because  effective  damping  of  fluid  film  bearings 
tends  toward  zero  for  half -frequency  oscillations,  the  rule  of  thumb  has  arisen 
that  the  threshold  of  instability  is  reached  at  about  twice  the  first  critical 
speed  of  the  rotor.  The  basis  for  this  rule  of  thumb  will  be  examined  in  more 
detail  in  subsection  3  where  there  is  discussed  a  procedure  for  using  critical 
speed  maps  for  predicting  threshold  of  whirl  instability. 

Having  briefly  discussed  the  qualitative  nature  of  whirl  instability  and  its  dis¬ 
tinction  from  critical  speeds,  we  will  now  proceed  to  show  the  nature  of  the 
linearized  analysis  required  to  determine  the  whirl  threshold  speed.  In  the 
illustrative  treatment  below,  we  shall  consider  only  the  translatory  mode  of  a 
symmetric  rotor  supported  on  two  identical  bearings.  The  more  general  case  of 
the  translatory,  conical  and  bending  modes  of  a  non-symmetric  rotor  will  be  con¬ 
sidered  in  subsection  2.  For  this  general  case,  determination  of  whirl  thresh¬ 
old  speed  requires  use  of  the  computer  program  developed  in  this  study.  For  the 
simpler  case  presented  below,  whirl  threshold  speed  can  be  determined  analytically. 

For  the  translatory  mode  of  a  rigid  symmetric  rotor  the  gravity  and  inertia 
(D'Alembert)  forces  are  equally  borne  by  the  two  bearings.  The  kinematic  rela¬ 
tionships  between  the  rotor  and  stator  bearing  surfaces  at  one  bearing  are  at 
every  instant  identical  as  those  of  the  other.  It  is  therefore  only  necessary  to 
consider  the  motion  of  one  journal  which  has  a  mass  equal  to  one-half  the  rotor 
mass.  Let  the  rotor  mass  be  2M  and  the  journal  center  amplitudes  be  x  and  y. 


At  any  given  rotor  speed  and  with  a  known  static  load  on  the  bearing,  the  journal 
center  occupies  a  certain  unique  equilibrium  position  relative  to  the  bearing 
center.  When  the  journal  whirls  around  its  equilibrium  position  in  a  small  orbit, 
hydrodynamic  bearing  forces  are  generated  in  the  bearing  fluid  film.  These  dy¬ 
namic  forces  can  be  expressed  in  a  linearized  form  by  expanding  the  film  forces 
into  a  first  order  Taylor  series.  With  the  bearing  fluid  film  dynamic  forces 
represented  by  and  F  and  with  the  dynamic  external  forces  on  the  rotor  repre¬ 
sented  by  W  and  W  as  shown  in  Fig.  1,  the  linearized  equation  of  motion  become: 
x  y 
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Fig.  1  Coordinate  System  for  Farces  and  Displacements 
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(1) 


where  x  and  y  are  the  whirl  amplitudes  measured  from  the  static  equilibrium  posi¬ 
tion,  t  is  time,  and  the  four  radial  stiffness  coefficients  and  the  four  radial 
damping  coefficients  are  computed  from  the  lubrication  equation  (Reynolds  equation) 
as  described  in  Ref.  2.  For  a  given  bearing  geometry  and  known  lubricant  proper¬ 
ties,  the  eight  coefficients  are  functions  of  the  bearing  load  and  the  rotor  speed 
and,  if  the  lubricant  is  compressible  like  a  gas,  they  are  also  functions  of  the 
whirl  frequency.  At  the  threshold  of  Instability,  x  and  y  are  pure  harmonic  mo¬ 
tions  and  can  be  conveniently  expressed  in  terms  of  the  whirl  frequency  vt 


x  «  x  cos  (vt)  -  x  sin  (vt) 
c  s 


y  ■  y  cos  (vt)  -  y  sin  (vt) 
c  s 


These  equations  can  also  be  expressed  in  complex  form: 

x  -  Re  I (x  +  ix  )  ei'l,t:'l  -  x  cos  (vt)  -  x  sin  (vt) 

»  c  s  J  c  s 

O) 

y  -  Re  {(y  +  iy  )  elvtj-  -  y  cos  (vt)  -  y  sin  (vt) 

where  Re  {  ^  indicates  that  only  the  real  part  of  the  complex  equation  is  applic¬ 
able.  For  convenience  the  Re  j  notation  is  dropped  and  Eq.  (3)  is  expressed: 

x  *»  (x  +  ix  )  eivt 
c  s 


y  =  (yc  +  tys)  e 


When  these  equations  are  used  in  the  analysis  their  complete  meaning  is  defined 
by  Eq.  (3).  Time  derivatives  can  also  be  carried  out  in  the  complex  notation, 
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bearing  in  mind  that  only  the  real  part  is  of  significance;  e.g. 


1 

at 


Re 


{<*« +  lx.>  *lvt} 


-  Re 


-  Re 


[(xc  +  ^  £  (eiVt)} 
<*c  +  iX8>  ^ 


■  -  vx  sin  vt  -  vx  cos  vt 
c  s 


By  differentiating  Eqs .  (4)  and  substituting  into  Eqs .  (1)  the  equations  of  motion 
can  be  written: 


M 

X 

■  -  z 

x  -  Z 

y  +  w. 

XX 

xy 

M 

y 

-  -  Z 

x  -  Z 

y  +  W, 

y* 

yy 

where : 


Z 


K  +  iv  B 

XX  XX 


xx 


-f  iyuo  B 


xx 


(Similarly  for  Z  ,  Z  ,  Z  ) 
'  xy  yx  yy 


(5) 


(6) 


V  *  v/co  (7) 

Here,  u)  is  the  angular  speed  of  rotation  and  Y  is  the  ratio  of  the  whirl  frequency 
to  the  rotational  frequency.  In  this  form,  the  equations  are  equally  valid  for  an 
Incompressible  and  a  compressible  lubricant.  In  matrix  form  Eqs.  (5)  can  be  ex¬ 
pressed  : 
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For  the  rotor-bearing  represented  by  Eq.  (8)  to  be  unstable,  it  is  necessary  and 

sufficient  that  Eq .  (3)  yield  a  non- zero  solution  for  the  amplitudes  x  and  y  when 

the  external  forces  W  and  W  are  zero,  i.e.  for  the  system  to  be  unstable  a  non- 

x  y 

trivial  solution  to  Eq.  (9)  below  must  be  found. 


(9) 


For  a  non-trivial  solution  to  exist,  the  determinant  of  the  matrix  must  vanish. 
Setting  both  the  real  and  imaginary  parts  of  the  determinant  to  zero  gives; 


A  •  A  +  iA  -  (Z  -  Mv2)(Z  -  Mv2)  -  Z__  Z  -  0 
c  s  xx  yy  j« 


(10) 


A  -  Re  fA]  ■  (K  -  Mv2)(K  -  Mv2)  -  K  K 
c  1  1  '  XX  yy  xy  yx 


V  (cuB  a*  -  U)B  cuB  )  *  0 
xx  yy  xy  yx 


(ID 


A  *  I  {A}  “  Y  f  (K  -  Mv2)  CUB  +  (K  -  Mv2)  ojB 
s  m  ^  J  l  xx  yy  yy  * 


K  cuB  -  K  cuB  ! 
xy  yx  yx  xyj 


(12) 
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Thus,  when  both  the  real  and  Imaginary  parts  of  the  determinant  vanish  simultan¬ 
eously  harmonic  motion  or  whirl  is  present  and  this  whirl  is  referred  to  as  rotor 
instability.  Every  combination  of  y  and  <0  satisfying  either  Eq.  (11)  or  Eq.  (12) 
can  be  represented  by  a  point  in  a  y -<d  plot.  Typically  the  locil  of  these  points 
appear  as  shown  in  Fig.  2. 

Eqs.  (11)  and  (12)  contain  two  unknowns,  the  whirl  frequency  ratio,  y,  and  the 

angular  speed  of  rotation,  CD.  For  incompressible  fluids,  the  eight  dynamic  film 

coefficients  are  Independent  of  y,  thus  Eqs.  (11)  and  (12)  may  he  solved  yielding 

2  2 

two  expressions  for  Mv  and  y  . 


Mv 


K  cuB  +  K  cuB  -  K  cuB  -  K  cuB 
xx  yy  yy  xx  xy  vx  vx  xy 

cub  +  cuB 

xx  yy 


(13) 


(K 


xx 


Mv2)(K  -  Mv2) 

In _ 


K  K 

xy  yx 


cuB  cuB  -  cuB  cuB 

xx  yy  xy  yx 


(14) 


2 

At  a  given  rotational  speed  and  load  and  with  the  eight  coefficients  known,  Mv 

may  be  calculated  from  Eq.  (13)  and  when  substituted  into  Eq .  (14)  gi-  s  the 

2 

whirl  frequency  ratio  squared  value  y  .  Two  methods  of  establishing  an  insta¬ 
bility  criteria  may  be  obtained  from  these  two  values.  First,  the  actual  mass 

2  2 

at  the  journal  bearing  may  be  substituted  into  the  Mv  term  to  solve  for  the  v 

2 

value  which  when  substituted  into  the  y  value  yields  the  apparent  threshold 

frequency.  If  this  threshold  frequency  is  identical  to  the  rotor  speed  at 

which  the  coefficients  were  based  then  the  threshold  speed  has  been  obtained. 

2 

The  second  method  substitutes  the  rotor  speed,  cu,  in  the  y  term  to  find  the 

7 

whirl  frequency  v  which,  when  substituted  into  the  Mv  term,  yields  a  critical 

mass  M£  for  onset  of  instability.  If  the  actual  journal  bearing  mass  M  is  less 

than  M  ,  then  the  system  is  stable;  if  M  is  greater  than  M  ,  the  system  is 
c  c 
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unstable.  The  mathematical  derivation  of  this  stability  criterion  is  given  in 
Kef.  4. 


From  compressible  fluids,  the  eight  radial  dynamic  fluid  film  coefficients  are 

functions  of  both  y  and  o>,  making  a  closed  form  solution  to  the  simultaneous 

Eqs.  (11)  and  (12)  impossible  and  Che  solution  is  most  conveniently  obtained 

graphically.  For  any  fixed  value  of  o>,  A  and  A  can  be  plotted  as  functions 

of  y  to  find  their  zero  points.  With  y  >  0  it  is  seen  that  A  has  one  zero 

point  and  A  has  up  to  two  zero  points  (only  true  in  this  simple  case).  The 
c 

calculation  is  repeated  for  several  values  of  au  and  the  results  for  the  vari¬ 
ous  zero  points  may  be  plotted  as  shown: 


Fig.  2  Loci  of  Roots  for  Real  and  Imaginary  Parts  of  Eq.  (10) 


The  intersection  of  the  two  curves  define  the  speed  at  which  instability  sets  in 
or  the  point  at  which  both  the  real  and  imaginary  parts  of  the  simplified  determ¬ 
inant  vanish  simultaneously. 
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Analysis  of  Stability  of  Arbitrary.  Non -Svante trie a I  Rotor 


Consider  the  arbitrary  rotor  shown  In  Fig.  3.  Assuming  that  the  two  ends  of  the 

rotor  are  free,  the  bending  moments  M*  and  M'  and  shear  forces  V'  and  V'  at 

atm  ym  xm  ym 

the  rotor  end  will  be  linearly  proportional  to  the  displacements  x^  and  y^  And  bend¬ 
ing  angles  0^  and  ©^  of  the  shaft  at  station  1  provided  these  displacements  and 
bending  angles  are  small  enough  such  that  linearized  analysis  applies.  Mathe¬ 
matically  we  can  write  this  linear  relationship  as 


(15) 


where  the  terms  a^,  ai2>etc',are  the  dynamic  Influence  coefficients  which  relate 
the  forces  at  station  m  to  the  motions  at  station  1.  These  coefficients  are 
functions  of  the  rotor  inertia,  rotor  flexibility,  and  the  dynamic  stiffness  and 
damping  coefficients  of  the  bearings  supporting  the  rotor.  Since  a^,  a12*etc-> 
are  dynamic  coefficients,  they  explicitly  depend  upon  the  rotor  speed  co  ar.d  the 
frequency  ratio  V  ■  v/a)  where  v  is  the  frequency  of  the  shaft  motions  being  con¬ 
sidered.  For  gas  bearing,  the  stiffness  and  damping  coefficients  contained  in 
all’  ai2,etc',are  als0  functions  of  oa  and  Y. 

For  symmetrical,  rigid  rotors  for  which  only  translatory  motions  x^  and  y^  are 
considered,  Eq.  (15)  above  reduces  to 
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X 


i 


STATION  n 


STATION  (n+ 1) 


Sign  Convention  for  Amplitude,  Slope,  8ending  Moment  ond  Shear  Force  Illustrated 
in  x  -  z  Plane 


Fig.  3 
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Thia  is  Che  situation  considered  in  the  previous  section,  and  we  see  by  comparisi 
of  Eq.  (15)  with  Eq.  (8)  that  the  coefficients  a^,  a^»  •jl*  *n^  *22  *r#  *lven  * 
the  analytical  expressions  j 

*11  " 

*12  “ 

*21  " 

*22 

For  the  more  general  case  of  an  arbitrary  flexible  rotor  supported  on  dissim¬ 
ilar  bearings,  the  coefficients  a^,  a^, etc.,  are  very  complex  functions  and 
usually  cannot  be  expressed  in  analytical  form  but  must  be  calculated  by  means 
of  a  computer  program.  Therefore,  one  of  the  principal  tasks  to  be  performed  by 
a  computer  program  written  for  analyzing  the  stability  of  arbitrary  rotors  is 
that  of  calculating  all  of  the  dynamic  influence  coefficients  as  a  function  of 
rotor  geometry,  rotor  speed,  frequency  ratio  Y  and  the  dynamic  bearing  coeffici¬ 
ents  supplied  as  input  to  the  program.  The  procedure  by  which  these  Influence 
coefficients  are  calculated  is  described  in  Ref.  3,  Appendix  VIII. 


(Zxx  ’  Mv2) 


Z 

xy 

(17) 

Z 

y* 

(Z  -  Mv2) 

yy 


Once  the  dynamic  influence  coefficients  are  calculated  for  an  appropriate  range 
of  a)  and  Y,  the  various  thresholds  of  stability  of  the  rotor-bearing  system  are 
determined  by  essentially  the  same  procedure  as  for  the  simple,  symmetrical  roto 
bearing  system  discussed  in  the  previous  section.  That  is,  thresholds  of  insta¬ 
bility  are  found  by  determining  the  values  of  co  and  Y  at  which  both  the  real  and 
the  imaginary  parts  of  the  determinant  of  the  matrix  of  dynamic  influence  coeffi 
cients  simultaneously  go  to  zero.  Again,  this  is  accomplished  by  plotting  on  a 
graph  of  Y  versus  a>  separate  curves  for  the  conditions  where  the  real  and  the 
imaginary  parts  of  the  determinant  go  to  zero,  and  determining  the  thresholds  of 
instability  from  the  intersections  of  these  curves. 


In  the  case  of  an  arbitrary  rotor  which  is  flexible  and  unsymmetric,  determining 
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the  real  end  imaginary  parts  of  the  determinant  of  the  matrix  of  influence  coef¬ 
ficients  cannot,  in  general,  be  accomplished  analytically  but  must  be  accomplished 
numerically  by  the  computer  program  for  analyzing  stability.  The  procedure  is 
essentially  one  of  having  the  computer  program  calculate  the  values  of  the  real 
and  imaginary  parts  of  the  determinant  over  a  specified  range  of  whirl  frequency 
ratios  at  a  fixed  rotor  speed.  The  program  then  determines  the  zero  points  of 
the  real  and  imaginary  parts  of  the  determinant  by  quadratic  interpolation  each 
time  a  change  in  the  sign  of  these  functions  occurs.  It  should  be  clearly  noted 
here  that  the  computer  program  itself  does  not  determine  thresholds  of  stability; 
this  remains  to  be  accomplished  by  the  designer  by  plotting  the  curves  of  the 
loci  or  the  real  and  imaginary  roots  of  the  stability  determinant. 

The  discussion  presented  above  serves  to.  describe  in  a  qualitative  way  the  nature 
of  the  process  of  determining  the  stability  of  an  arbitrary  rotor  and  how  this 
process  is  implemented  by  means  of  a  computer  program.  A  detailed  description  of 
the  computer  program  developed  in  the  present  project  for  the  purpose  of  analyzing 
rotor-bearing  stability  is  given  in  Appendix  I. 

1..  Determination  of  Rotor-Bearing  Stability  by  Means  of  Critical  Speed  Maps. 

For  a  flexible  rotor  with  two  translatory  and  two  angular  degrees  of  freedom, i .e 
*1'  ^i*  ®i  and  ®i  tlie  determinant  of  the  matrix  of  influence  coefficients  is  a 
higher  order  polynomial  in  V  and  ©  than  is  the  simple  determinant  from  Eq.  (9). 
Consequently,  there  are  many  real  and  imaginary  roots  of  the  determinant  and  a 
number  of  different  thresholds  of  instability, i .e .,  a  plot  such  as  that  shown  in 
Fig.  2  will,  for  an  arbitrary  flexible  rotor,  contain  many  curves  for  A_  *  0  and 
*  0  and  many  curve  intersections.  Each  different  intersection  or  threshold  of 
stability  will  correspond  to  a  different  mode  of  instability  motion,  e  .  g .,  trar.s  la- 
tory  whirl,  conical  whirl, etc.  In  general,  only  the  mode  occurring  at  the  lowest 
value  of  ©  is  of  interest  since  this  defines  the  practical  operating  limits  of 
the  system.  However,  the  presence  of  a  number  of  modes  of  instability  makes  the 
determination  of  the  lowest  mode  a  very  complex  and  cosily  task,  even  with  the 
aid  of  the  computer  program.  Hence,  it  is  desirable  to  find  methods  of  approxim¬ 
ately  determining  thresholds  of  instability  so  that  the  exact  process  of  determinin 
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such  thresholds  by  mesns  of  the  stebllity  program  can  be  accomplished  with  a 
minimum  of  costly  searching. 

One  convenient  approach  for  approximately  locating  the  instabilities  of  most 
rotors  can  be  had  by  use  of  a  critical  speed  map  for  the  rotor  in  question.  A 
critical  speed  map  shows  the  various  natural  or  resonant  frequencies  of  a  rotor 
as  a  function  of  a  single  value  of  spring  stiffness  K  assigned  to  the  bearings 
supporting  the  shaft.  A  typical  critical  speed  map  is  shown  in  Fig.  10.  For 
complex  rotors,  such  maps  are  usually  obtained  by  means  of  a  computer  program 
written  specifically  for  this  purpose. 

As  was  discussed  earlier  in  subsection  1.,  whirl  instability  can  be  viewed  as 
a  condition  of  undamped  resonance.  In  particular  Lund  (Ref.  4)  has  shown  that 
if  this  view  point  is  taken,  then  at  the  threshold  of  stability  the  fluid  film 
bearing  can  be  represented  by  a  single  effective  spring  stiffness  coefficient 
K  and  damping  coefficient  B  given  by 


,  f(K  -  K  )  (YoB  -  YdB  )  +  -i(K  YoB  +  K  YnB  ) 

K  .  JL/jr  +  x  )  .  Has _ 3x- _ £*- . -  yy. _ 2  ..**  v*  v*  *y 

2l  xx  yy'  A 

(18) 

B  -  t(VujB  +  YnB  )  -  A  (19) 

l  ax  yy 

where  A  is  expressed  by  the  following  quartlc  equation: 


A4  ♦ 


•f  (X  -  X  )2  +  X  X  -  7  (Vu)B  -  VccB  )2  -  Yob  YoB  I  A2 
4  xx  yy  xy  yx  4  xx  yy  xy  yx 


] 


Z(Kxx‘  V(Ya'6xx  -  V^Byy)  +  2  <Kxy  ^Byx  +  Kyx  Yi)Bxy  ] 


(20) 
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Sine*  B  is  real,  Eq.  (20)  has  only  two  solutions  for  A  which  are  of  equal  magni¬ 
tude  but  of  opposite  sign.  In  order  for  the  effective  bearing  damping  to  become 
zero  (corresponding  to  the  undamped  resonance  criteria)  at  some  value  of  V,  then 

A  must  be  positive  from  Eq.  (20)  since  (yoB  +  )  is  always  positive.  There- 

xx  yy 

fore,  only  the  real,  positive  value  of  A  is  used  to  define  the  stiffness  coef¬ 
ficient  in  Eq.  (18) . 


Another  method  of  obtaining  the  effective  stiffness  and  damping  bearing  coeffic- 

2 

ients  at  the  threshold  of  stability  is  by  replacing  MV  by  Z  in  Eq.  (9)  and  then 
solving  for  Z.  From  Eq.  (9): 


(Z  -  Z) 
xx 


y* 


”1  H 

<Zyy-Z)J  LyJ 


(21) 


Since  two  solutions  of  Z  are  obtained,  the  following  notation  is  used: 


*1-  Z2  “ 


major 


minor 


1 


-  Z  r  +  4Z  Z  I 
xx  yy  xy  yx  J 


(22) 


where : 


xx 


K  ♦  iYcuB 

XX  XX 


(Similarly  for  Z  ,  Z  ,  Z  ) 
xy  yx  yy 


Z  .  -  K  .  +  iYcoB 

major  major  major 


(Similarly  for  Z  .  ) 

J  minor 


The  condition  tor  the  lowest  threshold  speed  is  that  yuB 


minor 


*Sninor  and  aminor  are  e9uiva^ent  to  the  two  expressions  shown  in 
(19). 


0.  Thus  the 
Eqs .  (18)  and 
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In  arriving  ac  a  Kninor  *n<1  *  ®minor  v**u#»  t*'e  n®gfltive  sign  in  front  of  the 
radical  in  Eq.  (22)  is  generally  used.  An  evaluation  of  the  value  of  the  rad¬ 
ical  la  somewhat  difficult  since  the  coefficients  of  the  bearing  change  con¬ 
stantly  with  increasing  eccentricity.  Various  bearing  types  and  coefficients 
have  been  evaluated  and  this  sign  convention  has  proven  correct.  Thi*  sign 
convention  is  also  consistent  with  the  analysis  by  Lund  (Ref.  4).  ' ' '  , 

The  importance  of  Km£ncr  lies  in  the  convenience  with  which  it  can  be  used  in 
conjunction  with  a  rotor  critical  speed  map  to  determine  the  onset  of  insta¬ 
bility.  By  cross-plotting  the  minor  stiffness  of  the  bearing  as  a  function 
of  rotor  speed,  several  critical  speed  curves  may  be  crossed  before  the  rotor 
design  speed  is  reached.  The  threshold  speed  is  obtained  when  the  first  crit¬ 
ical  speed  corresponding  to  KB£nor>  divided  by  the  critical  whirl  frequency 
ratio  v,  gives  the  rotor  speed  at  which  the  Km^nor  and  y  were  evaluated.  The 
critical  whirl  frequency  ratio  is  defined  as  the  value  of  v/u)  at  which  B  eval¬ 
uated  from  Equation  (19)  goes  to  zero.  In  most  Instances,  the  critical  fre¬ 
quency  ratio  is  nearly  1/2,  and  since  ^m£nor  tends  not  to  change  very  rapidly 
with  rotor  speed,  this  gives  rise  to  the  rule  of  thumb  that  the  whirl  threshold 
speed  is  approximately  twice  the  first  critical  speed  of  the  rotor  bearing  system. 
The  critical  whirl  frequency  ratio  can  be  determined  from  Eq.  (14). 

The  specific  steps  Involved  in  determining  threshold  speed  from  a  critical  speed 
map  are  as  follows: 

(1)  Obtain  a  critical  speed  map  of  the  rotor  as  a  function  of  support 
stiffness. 

(2)  Calculate  the  eight  radial  bearing  coefficients  over  some  specified 
speed  range  which  one  estimates  will  contain  the  whirl  threshold  speed. 

(3)  Compute  the  and  y  values  from  Eqs .  (18)  or  (22)  and  (14),  for 

various  shaft  speeds. 
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(4)  For  each  of  Che  K  .  values,  enter  the  critical  speed  map  at  each 

nxnor 

valua  and  obtain  the  first  natural  frequency  or  critical  speed. 

(5)  For  each  shaft  speed  selected  divide  the  critical  speed  by  the  com¬ 
puted  whirl  frequency  ratio.  When  the  resulting  speed,  which  we  shall 
refer  to  as  apparent  threshold  speed,  equals  the  shaft  speed  at  which 
the  Km^nor  value  was  evaluated,  then  the  actual  threshold  speed  has 
been  determined. 

For  incompressible  lubricated  bearings  the  stiffness  and  damping  coefficients 
are  independent  of  the  whirl  frequency  ratio,  whereas  the  compressible  lubri¬ 
cated  bearing  coefficients  must  be  obtained  at  the  critical  whirl  frequency 
ratio.  In  general  these  gas  bearing  coefficients  may  be  obtained  at  v/o>  -  0.50 
with  sufficient  accuracy  since  the  critical  whirl  frequency  ratio  is  usually 
''lose  to  this  value. 

It  should  be  noted  that  use  of  a  critical  speed  map  together  with  values  of 
*minor  to  determine  the  stability  threshold  is  valid  for  rotors  of  ' trary 
shape  and  flexibility  but  does  rely  upon  the  bearings  supporting  the  rotor 
being  quite  similar,  since  the  critical  speed  map  is  plotted  in  terms  of  only 
one  value  of  bearing  stiffness.  However,  if  the  bearings  supporting  the  rotor 
are  not  too  dissimilar,  the  present  approach  will  serve  to  indicate  at  least 
approximately  the  critical  speed  and  critical  speed  ratio  by  using  values  of 
^minor  ^or  an7  one  t*'e  bearings.  A  more  exact  value  for  critical  speed  can 
then  be  determined  by  using  the  computer  program  developed  for  calculating 
rotor  bearing  stability. 

An  example  calculation  of  how  to  determine  whirl  threshold  speed  from  a  crit¬ 
ical  speed  map  is  given  in  the  next  section. 
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SECTION  III 


EFFECT  OF  THRUST  BEARING  STIFFNESS  AND  DAMPING  ON  ROTOR-BEARING  STABILITY 


1.  Discussion 

Since  a  thrust  bearing  exhibits  no  radial  stiffness  or  radial  damping  properties, 
the  threshold  speed  of  a  rotor  which  whirls  in  the  lateral  or  radial  mode  will 
not  be  affected.  However,  if  the  rotor  motion  is  angular  or  conical,  considerable 
restraint  can  be  imposed  on  the  rotor  by  the  angular  thrust  bearing  properties. 
Thus,  if  the  lowest  critical  speed  of  the  rotor  is  conical,  the  threshold  speed 
of  the  rotor  could  be  significantly  influenced  by  the  addition  of  the  thrust  bear¬ 
ing  coefficients.  The  computer  program  described  in  Appendix  I  of  this  report 
allows  for  the  inclusion  of  these  thrust  bearing  coefficients.  In  this  present 
section  we  will  consider  some  simplified  relations  which  enable  us  to  estimate 
the  effects  of  thrust  bearing  stiffness  on  the  conical  stability  of  a  symmetric 
rotor. 

First  let  us  define  the  axes  about  which  the  rotor  undergoes  angular  displacement 
while  whirling  in  a  conical  mode.  Consider  Fig.  4.  The  two  coordinates  which 
serve  to  define  the  angular  displacement  of  each  thrust  ber.ring  are  the  8  coor¬ 
dinate,  which  measures  rotation  about  the  y  axis  and  the  cp  coordinate  which  mea¬ 
sures  rotation  about  the  x  axis .  Both  of  these  coordinates  are  shown  about  point 
0  at  the  center  of  the  rotor  bearing  system.  These  angular  rotations  are  defined 
by  tan  9  -  dx/dz*  and  tan  cp  «  dy/dz*.  The  equations  defining  the  angular  dynamic 
coefficients  are: 


M 

-  G  9 

-  D 

d9 

-  G  cp  -  D 

dtp 

(23) 

X 

XX 

XX 

dt 

xy  xy 

dt 

M 

-  G  9 

-  D 

d0 

-  G  cp  -  D 

dg> 

(24) 

y 

xy 

yx 

dt 

yy  yy 

dt 

where  M  and  M  are  the  dynamic  moments  acting  on  the  rotor  resulting  from  the 
angular  displacements  and  velocities  9,cp  ,  d6/dt  and  dep/dt.  is  in  the  z 

plane  while  My  is  in  the  y-z  plane. 


*For  a  flexible  rotor,  local  values  of  dx/dz  and  dy/dz  serve  to  define  the  angular 
rotations  9  and  cp  at  local  points  along  the  rotor. 
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The  computer  program  written  to  analyze  rotor-bearing  stability  accepts  the  eight 

angular  stiffness  and  damping  coefficients  G^,  cdD^,  G  ,etc ..directly  as  bearing 

input  data  along  with  the  translatory  coefficients  K  ,  coB  ,  K  ,  etc.  However, 

xx  xx  xy 

in  order  to  gain  an  estimate  of  the  effect  of  thrust  bearing  angular  stiffness  on 
rotor-betring  stability,  it  is  convenient  to  try  to  use  the  critical  speed  map  ap¬ 
proach  for  determining  whirl-instability  that  was  described  earlier  in  the  report. 
To  use  this  approach,  it  is  necessary  to  relate  the  thrust  bearing  angular  coef¬ 
ficients  to  effective  translatory  stiffness  and  damping  coefficients  of  the  journal 
bearings  supporting  the  rotor.  For  a  symmetrical  rotor-bearing  system  this  can  be 


accomplished  by  the  following  geometrical  considerations.  Referring  to  Fig.  A  we 

see  that  the  restoring  moment  (M  ) ,  in  the  x-z  plane  exerted  by  each  journal  bear- 

^  J 

ing  in  response  to  a  static  angular  displacement  9  about  an  axis  through  the  point 


0  parallel  to  the  8  axis  is 


<»x>J 


-  K  L  0 

XX 


Similarly,  the  restoring  moment  (M^,  due  to  *ach  journal  bearing  is 


(M  ). 

y i 


*yxL  9 


On  the  other  hand,  the  restoring  moments  (M^)^  and  exerted  by  each  thrust 

bearing  due  to  the  angular  displacement  9  are  given  by  Eqs.  (23)  and  (2A),i.e., 


<Mx>I  ■  '“xx9  <27> 

<yT  ■  -  v  9  (28) 

Therefore,  we  see  that  for  angular  displacement  in  this  symmetrical  system,  the 

restoring  moments  of  thrust  bearings  can  be  represented  simply  by  adding  the 

2  2 

"effective"  stiffness  coefficients  G^/L  and  G^/L  respectively  to  the  existing 
coefficients  K^and  of  the  journal  bearings.  This  applies  as  well  to  all  of 
the  stiffness  and  damping  coefficients.  This,  for  angular  displacements  about 
axes  through  the  point  0,  the  restoring  moments  can  be  determined  by  assuming 
that  each  journal  bearing  has  effective  stiffness  and  damping  coefficients  (Kxx)e; 
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<BxxV  (Kxy)e'etc-*8lven  b* 


(K  ) 

m 

K 

+ 

G  /L2 

xx  e 

XX 

XX 

(K  ) 

m 

K 

+ 

G  /L2 

xy  e 

xy 

xy 

(K  ) 

m 

K 

+ 

G  /L2 

yxe 

y* 

yx 

(K  ) 

m 

K 

+ 

j  /L2 

yy  e 

yy 

yy 

(B  ) 

m 

B 

+ 

1)  /L2 

xx  e 

XX 

xx 

(B  ) 

at 

B 

+ 

D  /L2 

xy  e 

xy 

xy 

(B  ) 

m 

B 

+ 

D  /L2 

'  yx  e 

vx 

yx 

(B  ) 

m 

B 

+ 

D  /L2 

yy  e 

yy 

yy 

The  approach  of  taking  account  of  thrust  bearing  effects  on  conical  motions  of 
symmetrical  systems  by  means  of  adding  effective  stiffness  and  damping  to  existing 
journal  bearings  permits  one  to  use  the  critical  speed  map  method  for  estimating 
whirl  instability  as  described  in  the  previous  section.  This  would  be  accomplished 
by  using  the  total  effective  values  of  stiffness  for  conical  motion  of  the  rotor, 
as  defined  in  Eqs.  (29)  to  evaluate  ^  and  overall  damping  B  as  defined  by 

Eqs.  (18),  (19)  and  (2>>) .  Note  that  this  procedure  is  valid  for  conical  whirl 
instability  only.  Thrust  bearings  would  have  no  effect  on  the  translatory  whirl 
mode  and  their  angular  coefficients  should  not  be  included  in  the  evaluation  of 
K^inor  and  8  if  this  mode  of  instability  is  being  determined. 


Some  sample  calculations  were  performed  by  the  critical  speed  map  methoo  to  ex¬ 
amine  the  effect  of  independently  adding  principal  angular  stiffness,  G  and  G 

xx  yy 

and  principal  angular  damping  D  and  D  on  the  conical  stability  of  a  symmetrical 

xx  yy 

rotor.  The  conclusions  reached  were  that: 
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(1)  Addition  of  principal  angular  stiffnesses  G xx  and  G  does  not  change 
the  critical  value  of  frequency  ratio  for  instability  but  can  in¬ 
crease  the  threshold  speed  ai  for  instability  by  increasing  the  natural 
frequency  of  the  conical  mode. 

(2)  Addition  of  principal  angular  damping,  D  and  D  ,  does  not  change  the 

XX  yy 

natural  frequency  of  the  rotor-bearing  system,  but  does  decrease  the 
critical  frequency  ratio  for  instability  (i.e.,  the  frequency  ratio  at 
which  damping  B  goes  to  zero)  thereby  increasing  the  threshold  speed  co 
at  which  instability  occurs. 

In  the  sub-section  below  are  discussed  various  specific  design  examples  for  which 
the  effect  of  thrust  bearings  on  rotor-bearing  stability  are  investigated. 

2.  Sample  Calculations 

The  basic  tool  used  in  computing  the  threshold  speed  is  computer  program  PN400. 
With  this  program  it  was  possible  to  obtain  the  threshold  speeds  exactly  with 
and  without  the  thrust  bearing  effects.  However,  since  much  time  searching  for 
threshold  speeds  is  generally  required  when  using  this  program  alone,  the  shorter 
critical  speed  map  technique  for  locating  threshold  speeds,  outlined  in  this  re¬ 
port,  was  first  used  to  obtain  an  approximate  value  of  threshold  speed  with  an 
exact  solution  being  then  obtained  using  the  computer  program. 

Two  rotors  will  be  analyzed  to  determine  their  threshold  speeds  with  and  without 
the  addition  of  a  thrust  bearing  system.  These  two  rotor  models  are  shown  in 
Figures  5  and  6. 

The  first  rotor,  Figure  5,  weighs  approximately  7.2  lbs.  and  is  supported  by  two 
gas- lubricated  hydrodynamic  journal  bearings.  The  load  is  identical  on  each  bear- 
ing.  Air  at  100  F  and  ambient  pressure  (14.7  psia)  is  supplied  to  each  bearing. 
Thrust  bearing  surfaces  are  available  at  either  end  of  the  rotor. 

The  second  rotor,  Figure  6,  weighs  approximately  18.75  lbs.  and  is  supported  by 
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two  got- lubricated  hybrid  journal  bearings.  The  loads  are  different  for  each 
bearing.  A  thrust  bearing  is  shown  at  the  left  head  end  of  the  rotor.  This 
particular  model  was  designed  as  a  high  pressure  (BP)  spool  for  an  actual  high 
temperature  turbo-compressor  but  the  design  was  rejected  based  on  its  threshold 
speed.  No  consideration  was  given  to  the  thrust  bearing  effects  on  the  rotors' 
stability  behavior.  In  this  example  the  thrust  bearing  effects  will  be  computed. 

Tables  X  and  II  list  the  nondlmensional  dynamic  coefficients  for  the  hydrodynamic 
and  hybrid  gas- lubricated  Journal  bearing  designs  respectively.  Table  I  lists 
the  eight  stiffness  and  damping  coefficients  as  a  function  of  eccentricity  (or 
load),  bearing  number  A  and  whirl  frequency  ratio  v/o  for  a  length-to-diameter 
ratio  of  1.00.  For  a  constant  load  on  a  bearing  (W  constant)  and  varying  speed 
the  coefficients  must  be  obtained  by  double  interpolation  on  A  end  load  in  Table 
I  since  A  changes  with  speed.  Although  data  are  shown  for  v/cn  ■  0.3,.  0.5  and 
1.00,  only  the  v/o  *0.5  data  was  used  in  the  first  example. 

The  data  for  rotor  No.  1  are  given  below 

Rotor  Model  #1 

Journal  Bearing  Design  (Hydrodynamic) 

D  -  1.50  in. 

L  -  1.50  in. 

C/R  -  .001  in/in. 

Pg  ■  14.7  psia 

W  ■  3.62  lbs.  (each  bearing) 

U  *  2.7  X  10  ^  lb. sec/in2  (air  at  100®F) 

Thrust  Bearing  Design  (Hydrostatic) 

R  *  1.00  in. 

o 

Rj_  -  0.333  in. 

Pg  =  14.7  psia 

P  -  73.5 

s 

C  =  .0011  in. 


Inherent  compensation 


TABLE  I 

FLA IN  CYLINDRICAL  BEARING 
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Dimensionless  Operating  Parameters  for  Journal  Bearings: 


A 


12nuN 

P  (  C  ' 

a 


12 n(2 ■  7  Y  10'9)(N) 
14.7  (.001)2 


-a 

*=  6.95  *  10  " (N)  where  N  =  speed,  rps 

-4 

=  1.16  X  10  (N 1 )  where  N'  =  speed,  rpm 


P  LD 

3 


3.62 

(14.?)<1.5)<1.5) 

*  0.11 


Values  of  journal  bearing  stiffness  and  damping  coefficients  were  obtained  from 
Table  I  for  various  pertinent  values  of  A  at  v/u>  =  0.5  for  W  =  0.11.  Corresponding 
to  these  values  of  journal  bearing  coefficients,  values  of  K  .  were  calculated 
for  cue  journal  bearings  from  Eq .  (22)  neglecting  thrust  bearing  effects. 

I  he  critical  speed  map  for  rotor  #1  is  given  in  Figure  7.  This  figure  shows  the 
tirst  two  critical  speed  curves  as  a  function  of  support  stiffness.  Also  shown 
i  a  c.irve  of  the  minor  stiffness  for  the  bearings  calculated  in  the  manner  des- 

c  r  l  be J  above  . 


Whirl  threshold  speed  is  obtained  from  Figure  7  in  the  following  way.  First,  if 

we  draw  m  appropriate  vertical  line  on  Fig.  7  it  will  intersect  both  the  first 

<  ri'u  j  1  |  curve  and  the  curve  of  K  .  vs.  speed.  Let  us  denote  the  values 

minor 

.  !  .  p-  ■.  .1  .,t  Uii  :-i*  intersections  as  and  iDj,  respectively.  Now,  corresponding  to 
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these  values  of  ,  the  natur.il  frequency  of  oscillation  and  the  rotor  speed 

at  which  Kminor  is  calculated,  there  will  be  some  value  of  whirl  frequency  ratio 
Y  »  v/a)  calculable  from  Eq.  (14).  If  this  calculated  value  of  Y  is  exactly  equal 
to  then  we  have  determined  the  whirl  threshold  condition  and  will  be  our 

whirl  threshold  speed  while  v  /a>  will  be  our  whirl  frequency  ratio.  If,  say, 
Vc^U5T  s^ou^  prove  to  be  less  than  our  calculated  value  of  Y,  then  we  must  try 
drawing  another  vertical  line  slightly  further  to  the  right  and  redetermine 
If  v^/j)  Is  less  than  the  value  of  Y  calculated  from  Eq.  (14)  then  our  new  ver¬ 
tical  line  should  be  drawn  slightly  further  to  the  left.  Usually  Y  ~  0.5  so  we 
would  first  try  drawing  a  vertical  line  such  that  v^/cu^  =  0.5. 


In  the  example  shown,  neglecting  thrust  bearing  effects,  the  whirl  threshold  speed 
is  found  to  be  *  3300  rpra  while  the  critical  whirl  frequency  ratio  is  Y  *  vcMj. 


0.457. 


It  should  be  noted  that  in  using  this  critical  speed  man  approach,  the  curve  of 

K  ,  was  calculated  for  the  condition  v/u)  =  0.5  rather  than  for  the  condition 
minor 

v/u)  =  0.457.  For  greater  accuracy,  one  could  recalculate  Kminor  at  the  predicted 
critical  whirl  ratio  cf  0.457  and  redetermine  the  whirl  threshold.  Usually,  how¬ 
ever,  this  would  result  in  only  a  very  small  change  in  the  value  of  threshold 
speed  and  is  not  necessary. 

Next  let  us  consider  how  we  would  determine  the  threshold  speed  of  rotor  No.  1  in 
a  more  precise  manner  using  the  computer  program  PN400 .  In  using  the  program,  we 
can  take  advantage  of  the  fact  that  we  have  already  determined  the  threshold  speed 
in  an  approximate  manner  by  the  critical  speed  map  approach  as  described  above. 
Therefore,  with  the  computer  program,  we  look  for  the  roots  of  the  real  and  imag¬ 
inary  parts  of  Eq .  (10)  in  the  vicinity  of  a)  =  3300  rpm  and  Y  =  0.457. 

I  he  results  obtained  from  the  computer  program  are  shown  in  Fig.  8.  The  straighter 
Line  represents  the  locus  of  roots  of  the  imaginary  part  of  Eq .  (10)  ,  i  .  e  represent 
root-,  of  Kci  .  (12), while  the  more  curved  line  represents  the  locus  of  roots  of  the 
i  •  1 1  put  oi  Eq.  (10),i.e  , represents  roots  of  Eq  .  (11).  The  intersection  of  these 
m- r:  presents  the  condition  at  the  threshold  of  whirl  instability.  As  can  be 
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seen,  thin  occurs  at  a  running  speed  of  approximately  3570  rpm  and  at  a  frequa 
ratio  of  V  «  0.460.  j 

For  comparison,  the  whirl  threshold  solution  obtained  by  the  approximate  crit^ 
speed  map  approach  is  also  shown  in  Fig.  8.  Agreement  with  the  more  precise  c 
puter  solution  is  quite  reasonable. 

We  can  now  proceed  to  investigate  what  change  in  the  calculated  value  of  whirl 
threshold  speed  for  our  example  rotor  No.  1  will  result  if  we  consider  the  thr 
bearing  angular  stiffness  and  damping  in  our  calculations.  The  thrust  bearing 
for  this  rotor  are  hydrostatic  with  dimensions  listed  on  page  25.  In  the  case 
hydrostatic  thrust  bearings,  stiffness  is  the  only  bearing  characteristic  one 
usually  needs  consider  in  determining  the  influence  of  the  thrust  bearing  on  c 
icat  mode  of  whirl  instability.  This  is  so  for  the  following  reasons: 

(1)  Damping,  being  the  result  of  hydrodynamic  forces  rather  than  hydrosta 
forces, is  a  very  much  weaker  force  in  externally  pressurized  bearings 
than  is  the  stiffness. 

(2)  The  "effective"  angular  damping  in  hydrostatic  thrust  bearings  tends 
zero  at  Y  *  0.5  much  the  same  as  does  the  effective  radial  damping  in 
plain  cylindrical  journal  bearings.  Hence,  for  whirl  frequency  ratic 
near  0.5,  hydrostatic  thrust  bearing  damping  will  be  quite  ineffecti\ 
in  increasing  the  threshold  speed  for  conical  whirl  instability. 

To  determine  the  effects  of  the  hydrostatic  thrust  bearing  stiffness  on  our  ca 
culated  value  of  whirl  threshold  speed,  we  will  first  do  a  rough  calculation  t 
the  critical  speed  map  approach  and  then  do  a  refined  calculation  using  comput 
program  PN400.  Data  for  the  angular  stiffness  of  hydrostatic  thrust  bearings 
given  in  Figs.  13  through  16  in  the  next  section  as  a  function  of  the  dimensic 
feeding  parameter  •  As  discussed  in  the  next  section,  this  data  is  for  stat 
or  steady  state  displacement  of  the  bearing.  However,  in  most  cases  this  is  t 
plieable  to  dynamic  displacement  of  the  bearing  if  the  frequency  v  of  dynamic 
oscillation  is  low  enough.  To  determine  if  steady  state  data  is  applicable,  v 
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calculate  the  dimensionless  squeeze  number  a 


14.7(.0011)2 


-  .605  X  10"3v 

At 

Nt  -  3570  rpm 
so  *  373  rad /sec 
v  ■  (v/o3)co 

■  187  rad/sec 

Therefore 

CT  *  .114 

From  Fig.  18,  we  see  that  squeeze  number  of  c  =  0.114  is  sufficiently  low  such 
that  steady  state  values  for  bearing  stiffness  apply  (see  discussion  of  Figs.  17 
and  18  in  the  next  section) . 

To  determine  the  angular  stiffness  G  from  Fig.  16,  we  must  determine  the  dimension¬ 
less  feeding  parameter  A  .  Our  approach  is  to  set  A  at  the  value  yielding  maximum 
stiffness  i.e.,A  *»  1.0  .(This  optimum  value  for  A  can  be  achieved  for  our  bearing 

s  j  s 

by  the  proper  choice  of  na  in  A  .)  The  maximum  value  of  dimensionless  angular 

_  3 

stiffness  G  is 
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'  tSS9k$£ 


1  +  5 


1  +  2/3  5 


.12 


2  , 

Now,  our  bearing  is  inherently  compensated,  so  6  *  a  / cd  is  sufficiently  large 
(d  *  2a  and  c  «  a)  that  the  factor 


-  II  ±  o 

(1  +  2/3  62) 


(See  References  7  and  8) 


as  6  ->  «•  becomes  approximately 


in 

-i  + 1 

. 2  3 


Solving  for  G  from  the  expression  given  on  Fig.  16  we  have 


[‘ 


j"1 


tt  (R  2  -  R  2)  R  2  (P  -  P  ) 
o  i  o  s  a 


1  +  2/3  6^ 


12  tt  L  (  i  ) 2  ■  (  3  3 ) 2 1  (l)2  (73-3  -  14.7) 
( 3/2 ) (  00110) 


1  13  x  10  in  lb/radian 


As  described  earlier,  this  angular  stiffness  may  be  converted  to  an  effective 
increase  m  journal  bearing  stiffness  AKg  which  can  be  added  to  the  existing 
journal  hearing  stiffness  for  purposes  of  determining  whirl  tnreshold  speed 
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i 


-  1.15  X  104  lb/in. 


For  a  hydrostatic  thrust  bearing,  the  cross-coupling  terms  are  zero  and  the  stiff¬ 
ness  in  the  two  principal  planes  are  identical.  Since  Kminor  *s  increased  directly 

by  increases  in  K  and  K  ,  the  equivalent  increase  in  stiffness  £K  may  be  added 
xx  yy  e 

directly  to  the  *minor  obtained  previously  and  the  resultant  curve  for  Km^nor  is 
plotted  in  Fig.  7.  Using  our  graphical  approach  for  determining  whirl  threshold 
speed,  i.e.,  must  equal  y  calculated  from  Eq.  (14),  we  find  that  the  threshold 

speed  with  thrust  bearings  turns  out  to  be  at  approximately  6900  rpm.  Bearing 
stiffness  and  damping  coefficients  were  therefore  obtained  for  the  hydrodynamic 
journal  bearing  at  this  value  with  W  «  0.11  using  the  v/cu  *  0.50  data.  The  journal 
bearing  coefficients  were  then  submitted  along  with  the  actual  angular  thrust  bear¬ 
ing  coefficients  into  the  computer  program.  The  results  are  shown  plotted  in  Fig. 

9  which  yields  a  threshold  speed  of  7170  rpm  at  v/u>  *  0.477. 


Rotor  model  No.  1  is  an  idealized  model  which  serves  as  an  example  calculation  of 
how  thrust  bearing  stiffness  may  significantly  improve  the  stability  of  the  rotor 
to  fractional  frequency  conical  whirl  instability,  and  how  this  improvement  may  be 
calculated.  In  rotor  model  No.  2  we  have  an  example  of  an  actual  prospective  de¬ 
sign  for  a  high  temperature  turbocompressor  which  was  rejected  on  the  basis  of  its 
poor  stability  characteristics.  A  simplified  drawing  of  the  rotor  is  shown  in 
Fig.  6  and  a  brief  description  of  the  rotor  was  given  earlier  in  this  section.  Due 
to  the  overhung  nature  of  the  design,  the  lowest  mode  of  fractional  frequency 
whirl  instability  was  conical  in  nature.  However,  since  the  analytical  tools  were 
not  available  at  the  time  this  design  was  proposed,  no  account  was  taken  of  the 
effect  of  the  thrust  bearing  in  possibly  enhancing  the  rotor  stability.  We  shall 
examine  the  stabilizing  influence  of  this  thrust  bearing  in  what  follows  below. 


CONICAL  WHIRL  THRESHOLD  SPEED  WITH  THRUST 
BEARING  EFFECTS 
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Loci  of  Roots  for  Real  and  Imaginary  Parts  of  Stability  Determinant 


The  data  for  rotor  '-o.  2  ara  aa  follow-;- 1 


Thrust  Bearing  Design  (Hydrostatic) 

R  -  2.53  in. 

o 

R1  -  1.55  in. 

P  ■  18  psia 

a 

P  *  30.6  psia 

s 

T  -  1150°R 

8 


j 

i 

] 

J 

\ 

j 

j 

ft  -  2.472  x  105  in.2/sec2oR 
n  -  32  holes 

a  ■  .012  in.  (orifine  radius) 

Orlfic?  compensation 
Single  Plane  Emission 
(Ho.  2  bearing) 


M  ■  3.35  x  10  9  lb. sec. /in. 2 
C  -  .0015  in. 

n  -  40  holes 

a  ■  .030  in.  *  d/2  (inherently 

compensated) 

6  -  a2/cd  -  10 


Rotor  Model  #2 

Journal  Bearing  Design  (Hybrid) 

D  ■  3.00  in. 

1  ■  3.25  in. 

C  -  .0026  in. 

P  -  18  psia 

a 

P#  “  72  psia 

W  -  6.95  lb.  (No.  1  bearing),  11.80  lb. 

U  ■  3.40  x  10  9  lb.  sec. /in2 

T  -  760°R 


Table  II  gives  the  hybrid  gas-lubricated  jorrnal  bearing  data  for  an  L/D  *  1.08 
with  v/w  ■  0.50  at  Pfl/Pa  *  4.0.  The  hydrostatic  effect  parameter  Ag  for  the 

journal  bearing  design  with  orifice  compensation  (5  *  0)  and  single  plane  admis- 

* 

sion  is  given  by: 


A 


s 


6una2  Vrt 

P  C3  V  1  +  7 
8 


(6) (3.4  x  10~9) (32) ( .012) 2 


V 


2.472  x  10  x  760 


(72) (.0026)' 


-  1.025 
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The  data  in  Table  Hare  applicable  for  A  *1.0  with  the  only  variables  being  the 
hydrodynamic  effects  A  end  load  capacity  W.  Again  double  interpolation  is  required 
to  obtain  the  necessary  stiffness  and  damping  coefficients  at  a  constant  W  value. 
The  W  values  for  rotor  No.  2  are 


Journal  Bearing  Data 


W 


1 


P  LD 


a 


6.95 

18(3.25) (3) 


*  .0396  (bearing  J.1 


W2  P  LD 
a 


118 

18(3 .25) (3) 


1  .0672  (bearing  2) 

The  relationship  between  bearing  number  A  and  rotational  speed  N  is: 


-  2  38  n  10  ^  N  (where  N  =  rps) 

■  396  >  10  '  N'  (where  N'  =  rpm) 
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The  beering  data  shown  in  Table  III  are  obtained  from  Table  II,  by  interpolation, 
for  “  .0396  and  “  .0672.  The  values  of  Kninor  given  in  Table  III  are  calcu¬ 
lated  from  Eq.  (18)  or  Eq.  (22).  The  values  of  y  •*  v/<o,  the  frequency  ratio  at 
which  effective  damping  gees  to  zero,  are  calculated  from  Eq.  (14). 

The  critical  speed  curves  for  rotor  No.  2  are  plotted  in  Figure  10.  Also  plotted 

is  the  curve  of  K  ,  vs.  speed.  In  this  case  K  .  varies  hardly  at  all  with 
minor  r  minor 

speed,  which  makes  calculation  of  the  whirl  threshold  speed  cu^,  very  easy.  Noting 
that  for  instability  vcMj  must  equal  V  (which  in  this  case  is  very  nearly  exactly 
0.5)  and  noting  that  vc  can  now  be  given  by  the  intersection  of  the  curve  of  Kminor 
with  the  lowest  critical  speed  curve,  we  have 

■  14,750  rpm  (from  Figure  10) 

■  v  /0.5  *  29,500  rpm 

To  obtain  a  more  exact  calculation  for  a^,  we  use  computer  program  PN400  with 
bearing  input  data  obtained  from  Table  III.  The  results  obtained  are  shown  in  Fig. 
11.  The  approximate  solution  obtained  above  from  the  critical  speed  map  approach 
is  shown  as  a  circle  in  Figure  11.  The  solution  curves  for  the  real  and  imaginary 
parts  of  the  stability  determinant  being  equal  to  zero  (the  solution  curves  for 
Eq.  (11)  and  (12))  are  very  nearly  parallel  for  this  example,  making  exact  deter¬ 
mination  of  the  whirl  threshold  speed  quite  difficult.  As  can  be  seen,  the  whirl 
threshold  speed  lies  somewhere  in  the  range  of  29,000  to  29,375  rpm  and  the  whirl 
frequency  ratio  is  between  0.499  and  0.504. 


Next  we  compute  the  effect  of  the  thrust  bearing  on  this  whirl  threshold  speed. 

Using  the  thrust  bearing  design  data,  the  squeeze  number  c  evaluated  at  the  pre 

vious  calculated  value  of  v  ,  the  oscillation  frequency  at  whirl  threshold  is 

c 


18  ( .0015)2 

=  6.0 
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HYBRID  GAS  BEARING  DATA 
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The  restrictor  coefficient  A  equals: 


A 


s 


6una2  V  &t 
P  C3  V 1  +  52 

3 


.  (6) (3.33  X  10~9)(40)(.03)2  V2.472  x  105  x  1150 

<30.6)(.0015)3  Vl  +  Id2 


-  11.80 

This  v*lue  of  Ag  is  significantly  higher  than  the  optimum  value  for  maximum  bear¬ 
ing  stiffness  (see  Figure  15) .  This  is  because  the  operating  condition  being 
analyzed  is  the  one  at  which  stability  of  the  rotor  is  poorest  and  is  not  the 
operating  condition  for  which  the  bearing  was  designed.  The  curve  of  dynamic 
stiffness  shown  in  Fig.  17  pertain  to  a  value  of  Ag  =  2.5  which  is  optimum  for 
a  thrust  bearing  of  radius  ratio  of  R^/R^  ■  1.5.  From  this  curve  we  see  that 
such  a  bearing,  operating  at  a  squeeze  number  of  o  -  6.0  would  be  well  within 
the  range  where  steady  state  stiffness  data  apply;  i.e.,is  operating  at  an  os¬ 
cillation  frequency  which  is  well  below  the  frequency  at  which  squeeze  film 
effects  become  important.  Physically,  there  is  no  reason  to  expect  that  this 
situation  would  be  altered  significantly  when  the  bearing  is  operated  at  higher 
Ag  (i-e., lower  supply  pressure).  Therefore,  we  would  readily  expect  that  our 
example  bearing  at  Ag  «  11.8  and  a  =*  6.0,  would  be  in  the  "steady  state"  region 
where  Fig.  15  could  be  used  to  calculate  bearing  stiffness.  From  Fig.  14,  at 
P  /P  *  1.70,  we  determine  that 

S  3 


1  +  2/3  6 

G  =  (.07)  S2  .  (  Q7)  .LI.  2/3  (10) 2 

1  +  t  1  +  (10)2 
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.0469 


TT(R  2-  R.2)  R  2  (P  -  P  )  G 
o  i  o  s  a 


z  r(2.53)2  -  q.55)21  (2.53)2  (30.6  -  18H.0469) 

.0015 


=  3.02  y  104  In. lb/rad. 


Since  we  have  only  one  thrust  bearing,  we  must  appropriately  divide  its  effect 
among  the  two  journal  bearings  to  use  the  effective  stiffness  -  critical  speed 
map  approach  for  approximately  calculating  the  effect  of  the  thrust  bearing  on 
whirl  threshold  speed.  The  distance  of  journal  bearings  1  and  2  from,  the  c.g. 
of  the  rotor  are  =  4.72”  and  -  2.78"  respectively.  By  analogy  with  Eq . 
(29),  which  gives  the  relationships  for  adding  an  effective  stiffness  of  one 
thrust  bearing  to  one  journal  bearing,  we  can  infer  that  a  correct  approach  for 
■’dding  an  effective  stiffness  of  one  thrust  bearing  ot  angular  stiffness  G  to 
two  journal  bearings  of  spans  and  would  be  by  the  relationship 


Therefore,  the  effective  radial  stiffness  to  be  added  to  each  journal  bearing 
won  1  d  be 


K  = 
e 


Li  +L2 


3.02  y  10 
30 


---  1010  lb/in 
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This  value  is  seen  to  be  negligible  compared  with  the  values  of  stiffness 
associated  with  the  journal  bearings  themselves.  Hence,  in  this  case,  the 
thrust  bearing  would  have  negligible  effect  on  conical  motions  of  the  shaft. 

The  above  conclusion  reached  for  rotor  NoV  2  probably  pertains  to  most  rotor 
designs.  In  general,  such  rotors  will  be  designed  so  that  the  span  between 
journal  bearings  is  efficiently  great  that  these  bearings  will  exert  a  much 
greater  restraining  force  on  conical  motions  of  the  shaft  than  will  the  thrust 
bearings.  However,  this  need  not  always  be  true.  It  is  easy  to  conceive  of 
systems  where  the  thrust  loads  dominate  and  where  journal  bearings  are  needed 
only  to  locate  the  shaft.  In  such  systems  a  designer  could  take  advantage  of 
the  now  available  technology  for  including  thrust  bearing  effects  in  whirl 
threshold  speed  calculations,  and  design  the  rotor  such  that  the  thrust  bear¬ 
ing  could  assume  the  responsibility  of  restraining  conical  motions  of  the  rotor. 
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SECTION  IV 


ANGULAR  STIFFNESS  COEFFICIENTS  FOR  HYDROSTATIC  THRUST  BEARINGS 


Ttie  hydro-'iatic  bearing  data  presented  in  this  section  are  based  on  an  analysis 
by  Lund  described  in  Reference  [s].  Basically,  the  analysis  involves  the  assump¬ 
tion  that  discrete  feeding  holes  arranged  in  an  annular  thrust  bearing  (see  F;g. 
12)  can  be  represented  approximately  by  a  continuous  "line  source"  of  feeding  as 
if  t he  bearing  were  fed  by  a  continuous  circular  slot  rather  than  by  discrete 
holes.  l'his  approximation  tends  to  be  absolutely  correct  in  the  limit  as  the 
number  of  feeding  holes  approaches  an  infinite  number,  but  tends  to  somewhat  over¬ 
estimate  bearing  Load  and  stiffness  for  bearings  having  a  finite  number  of  holes. 

A  line  source  correction  factor  to  account  for  this  overestimation  has  been  worked 
out  by  comparing  a  detailed  discrete  feeding  analysis  with  the  simpler  line  source 
analysis  The  correction  factor  to  be  applied  depends  on  the  product  n?  and  d/?D 

where  n  is  the  number  of  feeding  holes,  ?  -  1/2  log^R^/R^),  d  is  the  diameter  of 

the  feeding  holes,  D  *  2  v  R  R  ,  R  and  R,  are  the  outer  and  inner  radii  of  the 

o  l  o  i 

thrust  bearings  having  a  practical  number  of  holes,  the  correction  factor  to  be 
applied  is  2/3, i.e.,  the  line  source  analysis  overestimates  the  stiffness  by  50%. 
This  correction  factor  has  been  built  into  the  design  curves  presented  In  this 
sec  tier.  .  Figure  L9  shows  a  curve  of  n?  vs.  d/?D  giving  the  recommended  size  and 
number  of  feedings  for  which  this  correction  applies.  If  more  feeding  holes  are 
used,  stiffness  will  be  higher  than  that  predicted  by  the  design  charts  in  this 
ret  ion  and  conversely. 

The  stiffness  of  compressible  hydrostatic  bearings  is  dependent  upon  the  frequency 
of  excitation  v  of  the  bearing  if  this  frequency  is  high  enough.  This  is  because, 
at  large  values  of  v,  the  gas  in  the  bearing  film  tends  to  be  dropped  ano  com¬ 
pressed  rattier  than  squeezed  out  of  the  film,  and  this  contributes  to  a  higher 
level  of  bearing  stiffness.  This  effect  is  illustrated  in  figures  17  and  18  where 

dimensionless  stiffness  G  =  CG/Tr(R  -  R.  )  R  (P  -  P  )]  is  plotted  vs .  the 

o  t  o  s  a 

dimensionless  excitation  frequency  a  (known  as  squeeze  number) 


12..  v 

P 


) 
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r.  holes 
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as  can  he  seen,  for  small  a,  the  stiffness  remains  constant  at  its  steady  state 
value  independent  of  excitation  frequency.  However,  once  v  becomes  large  enough, 
the  bearing  stiffness  begins  to  increase  with  v  due  to  the  above  mentioned  squeeze 
i i In  effect. 

In  general,  for  rotors  supported  on  hydrostatic  bearings,  the  running  speed  will 
be  in  the  range  such  that  an  excitation  frequency  v  =  tn/2  will  lead  to  squeeze 
numbers  sufficiently  low  such  that  steady  state  bearing  stiffness  data  will  apply. 
This  is  generally  so  because  if  shaft  speeds  are  high  enough  such  that  compress¬ 
ibility  becomes  significant  in  the  squeeze  film  effect,  then  the  shaft  speed  will 
tend  to  be  high  enough  such  that  hydrodynamic  bearings  could  be  used  for  load 
support.  Figures  17  and  18  will  be  used,  therefore,  primarily  to  define  the  lim¬ 
its  below  which  steady  state  bearing  data  apply  rather  than  being  used  as  a  source 
of  dynamic  stiffness  data. 

In  Figures  13  through  16  are  "resented  dimensionless  angular  stiffness  G  vs.  the 
feeding  parameter  A 

s 


P  C3  Vl  +  62 


s 


(31) 


tor  radius  ratios  of  R  /R.  =  1.25,  1.5,  20,  and  3.0.  The  dimensionless  stiffness 

O  1 

—  A  2  2  2 
G  is  p re -mu  It  ip l led  by  the  factor  1+6  /(l  +  2/3  6  )  where  5  =  a  cd ,  known  as  the 

2 

inherent  compensation  factor,  gives  the  ratio  of  orifice  area  rra  to  inherent  re¬ 
striction  area  rfcd  (see  Figure  12).  Curves  are  provided  foi  different  values  of 
pressure  ratio  P  /P  .  Usually,  the  hydrostatic  thrust  bearings  are  designed  for 

S  d 

a  value  of  A  which  provides  near  maximum  stiffness;  i.e.,  Ag  is  in  the  range  1.0 
to  4  0,  depending  on  radius  ratio.  Hence,  the  dynamic  curves  in  Figs.  17  and  18 
are  presented  for  near  optimum  values  of  A  . 


Onlv  stiffness  data  is  provided  for  the  hydrostatic  bearings.  This  is  for  the 
following  reasons: 
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(1)  Damping  in  hydrostatic  bearings  results,  essentially,  from  hydrodynamic 
forcea  which  are  usually  much  weaker'  than  the  hydrostatic  forces  produc¬ 
ing  bearing  stiffness. 

(2)  More  importantly,  in  hydrostatic  thrust  bearings,  ”ef fective"  bearing 
damping  tends  to  zero  for  wobbling  or  nutating  nodes  of  motion  which 
occur  at  half  the  rotational  speed  of  the  thrust  runner.  In  examining 
the  stability  of  rotors,  it  is  invariably  the  influence  of  precisely 
this  kind  of  thrust  bearing  motion  on  rotor  dynamics  that  we  are  at¬ 
tempting  to  calculate.  Hence  in  such  calculations,  it  is  quite  reason¬ 
able  to  presume  that  the  effects  of  thrust  bearing  damping  will  be  zero 
or  near  zero  and  only  consider  the  effects  of  thrust  bearing  stiffness. 


Fig.  13  Hydrostatic  Thrust  Bearing,  Rq/R.  =  1.25 

Dimensionless  Angular  Stiffness  vs.  Restrictor  Coefficient 
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Fig.  14  Hydrostatic  Thrust  Bearing,  Rq/R^  =  1.5 

Dimensionless  Angular  Stiffness  vs.  Restrictor  Coefficient 


Fig.  16 


i-J-l  1  i 


Stiffness  for  the  Hydrostatic  Thrust  Bearing 
Inherent  Compensation 
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Fig.  18  Dynamic  Angular  Stiffness  for  the  Hydrostatic  Thrust  Bearing 
Rq/R;  =  3.  Inherent  Compensation 


10"3  10'2  TO-1  1 


d/CD 


Fig.  19  Minimum  Number  of  Feeder  Holes 
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APPENDIX  I 


COMPUTER  PROGRAM 

THE  THRESHOLD  GF  INSTABILITY  OF  A  FLEXIBLE  ROTOR  IN  FLUID  FILM  BEARINGS 


This  section  describes  the  rotor  stability  computer  program  PN400:  "The  Threshold 
of  Instability  of  a  Flexible  Rotor  in  Fluid  Film  Bearings".  This  program  is  quite 
similar  to  the  computer  program  PN0017  given  in  Reference  1. 

The  new  program  described  in  this  volume  is  a  modification  and  extension  of  the 
previous  version.  The  amount  of  plotting  has  been  significantly  reduced  in  that 
a  quadratic  interpolation*  routine  now  automatically  performs  much  of  the  plotting 
previously  required  to  find  the  zero  point  solutions  for  both  the  real  and  imagin¬ 
ary  parts  of  the  complex  determinant.  This  program  also  now  accepts  the  angular 
thrust  bearing  properties  and  is  also  set  up  to  handle  both  liquid-lubricated  and 
gas-lubricated  bearings.  In  addition  a  more  refined  rotor  model  is  employed  which 
describes  the  rotor  in  terms  of  stiffness  and  mass  diameters,  lengths,  concen¬ 
trated  values  of  mass,  and  mass  moments  of  inertia.  The  new  rotor  model  is  dis¬ 
cussed  in  Reference  3 . 


Analysis 

The  rotor  analysis  is  described  in  detail  in  Reference  3,  Appendix  VIII.  Con¬ 
sider  the  rotor  shown  in  Fig.  3  and  assuming  the  two  ends  of  the  rotor  to  be 

free,  the  bending  moment  and  the  shear  force  at  one  rotor  end  are  directly  pro¬ 
portional  to  the  unknown  amplitude  and  slope  at  the  opposite  rotor  ^nd .  The 

proportionality  coefficients  are  the  dynamic  influence  coefficients  which  are  .  - 

v 

computed  by  the  program,  end  the  final  equation  is  given  in  the  .following  matrix 
form  (see  Eq .  (H.51)  in  Appendix  VIII  of  Reference  3  neglecting  magnetic  influence 
coefficient  matrix)  . 


Second-order  polynomial  interpolation.  See  Reference  6. 


61 


Here,  V'  and  V'  are  the  shear  force  components  at  one  rotor  end  (station  m) , 

’  xm  ym  r  '  ’ 

and  M'  ^  are  the  bending  moment  components  at  the  same  location,  and  x^,  y^, 

and  rp^  are  the  unknown  amplitudes  and  slopes  at  station  1.  The  a's  are  the 

dynamic  influence  coefficients  which  depend  on  the  whirl  frequency  ratio  and  the 

rotor  speed.  They  include  the  effect  of  rotor  inertia,  rotor  flexibility  and 

the  dynamic  bearing  coefficients.  Since  the  ends  of  the  rotor  are  free,  V'  ** 

V'  =  M'  =  M'  ®  0.  Thus,  x,,  y, ,  0,  and  m,  are  only  different  from  zero 

ym  xm  ym  1  1  Y1  1 

when  the  determinant  of  the  matrix  of  influence  coefficients  is  zero  which,  then, 
determines  the  threshold  of  instability: 


A  =  A  +  iA 

c  s 


'll 

a12 

a13 

a14 

'21 

a22 

a23 

a24 

'31 

a32 

a33 

a34 

\  1 

a42 

a43 

a44 

(33) 


For  each  given  angular  speed,  co,  of  the  rotor,  the  program  calculates  A^  and  A 
for  specified  values  of  the  whirl  frequency  ratio  y  =  v/'JU  and  determines  the 
zero  points  of  A^  and  Ag  within  the  given  range.  The  results  can  be  plotted 
directly  to  find  the  threshold  of  instability. 
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COMPUTER  INPUT 


Card  1  (72H) 

Any  descriptive  text  may  be  given  to  Identify  the  calculation. 

Card  2  (1415) 

This  is  the  "control  card"  whose  values  control  the  rest  of  the  input.  The  con¬ 
trol  card  has  10  items: 

1 .  NS,  the  number  of  rotor  stations,  see  "Rotor  Data"  (1  <  NS  <  100) 

2 .  NB,  the  number  of  bearings  (1  <  NB  <  10) 

3 .  NPST.  The  program  provides  for  including  the  response  characteristics  of 
the  bearing  pedestals.  Whan  it  is  desired  to  include  this  effect,  set  NPST  *  1 
and  give  the  required  input  data  as  explained  later.  When  NPST  ■  0,  the  pedes¬ 
tals  are  assumed  to  be  rigid  and  no  pedestal  data  are  required  in  the  input. 

4 ■  NANG.  The  dynamic  forces  of  the  fluid  film  in  the  journal  bearings  resist 
both  radial  motion  and  angular  motion.  In  most  cases,  only  the  radial  forces 
are  significant  in  which  case  NANG  is  set  equal  to  zero.  However,  long  jour¬ 
nal  bearings  and  especially  thrust  bearings  may  exert  considerable  constraint 
on  the  angular  motion  of  the  rotor.  When  it  is  desired  to  include  this  effect, 
set  NANG  =  1  and  give  the  required  input  data  as  explained  later. 

5 .  NFR  is  the  number  of  whirl  frequency  ratios  specified  in  the  input  list 

below,  see  explanation  later  (3  <  NFR  <  50). 

6 .  INC.  If  the  bearing  lubricant  is  compressible  (gas  bearings),  set  INC  =  1 
in  which  case  the  dynamic  bearing  coefficients  must  be  specified  for  each  whirl 
frequency  ratio  in  the  input.  If  INC  =  0,  the  lubricant  is  incompressible  and 
the  dynamic  bearing  coefficients  are  independent  of  whirl  frequency. 
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7 .  NIT.  As  discussed  previously,  the  program  searches  for  the  zero  points 
of  the  instability  determinant.  To  do  this,  the  instability  determinant  is 
computed  as  a  function  of  the  whirl  frequency  ratio.  The  input  gives  NFR 
values  of  the  whirl  frequency  ratio  in  sequence,  and  the  determinant  is  cal¬ 
culated  for  each  of  these  values.  In  addition,  each  interval  is  subdivided 
into  NIT  increments  and  the  determinant  is  also  computed  at  these  intermedi¬ 
ate  frequency  ratio  values.  Whenever  the  program  detects  a  change  in  sign 
of  the  determinant  between  two  consecutive  calculations,  it  uses  quadratic 
interpolation  to  find  the  accurate  zero  point. 

By  subdividing  the  intervals  in  the  given  frequency  ratio  list,  the  list  can 
be  shortened  without  loss  of  accuracy.  Furthermore,  when  the  bearing  lubri¬ 
cant  is  a  gas  it  is  necessary  to  provide  the  dynamic  coefficients  for  each 
specified  frequency  ratio.  Thus,  to  minimize  the  input  data  it  is  desirable 
to  keep  this  number  low.  Then  the  program  automatically  calculates  the  coef¬ 
ficients  for  the  intermediate  points  by  quadratic  interpolation. 

NIT  should  be  equal  to  or  greater  than  1.  When  NIT  «  1,  no  subdivision  takes 
place  . 

8  METR.  This  item  should  always  be  zero.  It  is  included  for  diagnostic 
purposes  which,  however,  is  of  no  value  to  the  general  user  of  the  program. 

9 .  NCAL,  specifies  the  number  of  rotor  speed  ranges.  For  each  speed  range 
it  is  necessary  to  give  input  data  for  the  dynamic  bearing  coefficients. 

There  is  no  limit  for  NCAL. 

10 •  IMP .  If  INF  =  0,  the  present  set  of  input  data  is  followed  by  a  new  set, 
starting  from  Card  1.  If  IN?  *  1,  this  is  the  last  set  of  input. 

Card  3  (1P6E12  .5) 

2 

l ■  YM  is  Youngs  modulus  E  for  the  shaft  material  in  Ibs/inch  .  If  E  actually 
changes  along  the  rotor  it  should  be  noted  that  the  program  only  uses  E  in  the 
product  El  where  1  is  the  cross-sectional  moment  of  inertia  of  the  shaft.  Sine 
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fV  - 

L  °  stiff  1  J 


where  (d  ) 


,  is  the  outer  9haft  diameter  and  d.  is 

o'  _  i 


stiff 


the  inner  diameter,  any  variation  in  E  can  be  accounted  for  by  changing  (dQ) 
(see  "Rotor  Data") . 


stiff 


2 ■  DENST  gives  the  weight  density  of  the  shaft  material  in  lbs/inch  .  The  pro¬ 
gram  converts  it  into  the  mass  density  o  *  DENST/386 .069 .  If  the  density  actually 
changes  along  the  rotor  it  should  be  noted  that  the  weight  of  the  shaft  per  unit 


length  is  o 


4  r<d 

4  o 

w 


y 


mass 


where  (d  ) 
o 


is  the  outer  shaft  diameter  (see 


mass 


"Rotor  Data").  Thus,  (d  ) 


can  be  changed  to  absorb  the  changes  in  density. 


mass 


3 ■  SHM  gives  the  product  oG  where  G  is  the  shear  modulus,  lbs/inch  ,  and  a  is 
the  shape  factor  for  shear  (for  circular  cross-sections:  a*  0.75). 


Rotor  Data  (8E9  .2) 


The  rotor  is  represented  by  a  number  of  stations  connected  by  shaft  sections 
of  uniform  diameter.  Thus,  rotor  stations  are  introduced  wherever  the  shaft 
diameter  changes  (or  changes  significantly).  Also,  there  must  be  a  rotor  sta¬ 
tion  at  each  end  of  the  rotor,  at  each  bearing  centerline  and  at  any  thrust 
bearing.  Furthermore,  a  rotor  station  is  introduced  wherever  the  shaft  has  a 
concentrated  mass  which  cannot  readily  be  represented  in  terms  of  an  inner  and 
outer  shaft  diameter  (impellers,  turbine  wheels,  alternator  poles,  and  so  on). 

In  this  way  the  rotor  is  assigned  a  total  of  NS  stations  (card  2,  item  1)  which 
are  numbered  consecutively  starting  from  one  end  of  the  rotor.  There  can  be  a 
maximum  of  100  stations.  Eaci:  station  can  be  assigned  a  concentrated  mass  m 

with  a  polar  mass  moment  of  inertia  I  and  a  transverse  mass  moment  of  inertia 

P 

(any  of  these  quantities  may,  of  course,  be  zero).  Also,  each  station  can 

be  assigned  a  shaft  section  with  which  it  is  connected  to  the  following  station. 

This  shaft  section  has  a  length  L,  and  outer  diameter  (d  )  ,  an  outer  diameter 

0  stiff 


(do> 

mass 
ify  the 
inertia 


and  an  inner  diameter  d,.  The  outer  diameter  (d  ) 

l  o 


is  used  to  spec- 


stiff 


stiffness  of  the  shaft  section  such  that  the  cross-sectional  moment  of 


of  the  shaft  is 


TT 
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(d  ) 
u 


-  d. 


stiff 


and  the  shear  area  is: 
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The  outer  diameter  (d  ) 

o 


mass 


mass  of  the  shaft  such  that  the  mass  per  unit  length  is: 
where  0  is  the  mass  density  (see  card  3,  item  2). 


is  used  in  calculating  the 


2l 


f  [«, 


-  d 


mass 
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In  the  computer  input  there  must  be  a  card  for  each  rotor  station  (NS  cards)  . 
Each  card  specifies  the  7  values  for  the  station: 


1  .  The  concentrated  mass:  m,  lbs.  (may  be  zero), 

2.  Tne  polar  mass  moment  of  inertia  of  the  station  mass:  I 

—  P 

zero)  . 

3 .  The  transverse  mass  moment  of  inertia  of  the  station  mass 
be  zero)  . 


lbs- inch  (may  be 


I  lbs- inch  (may 


4 .  The  length  of  the  shaft  section  to  the  next  station:  f,  inch  (may  be  :iero) . 
For  the  last  station,  set  l  =  0. 


5 ■  The  outer  diameter,  (d  )  of  the  shaft  section,  inch.  (d  )  is  used 

°  stiff  °  stiff 

in  calculating  the  stiffness  of  the  shaft  section,  (d^)  ¥  0.  For  the  last 

station,  set  (d  )  »  1.0.  stiff 

°  stiff 

b .  The  outer  diameter,  (d  )  of  the  shaft  section,  inch  (may  be  zero). 

mass 

(d  )  is  used  in  calculating  the  mass  of  the  shaft  section.  For  the  last 

mass 

station,  set  (d  )  =  0 . 

mass 


The  inner  diameter, 
both  in  calculating  the 
last  station ,  set  d ,  = 

l 


d,  of  the  shaft  section,  inch 

l 

stiffness  and  the  mass  of  the 


(may  be  zero) . 
shaft  section. 


d.  is  used 
l 

For  the 


bb 


0. 


Bearing  Stations  (1415) 


The  rotor  station  numbers  at  which  there  are  bearings,  are  listed  in  sequence. 
This  includes  both  journal  bearings  and  any  thrust  bearings.  There  car.  be  up 
to  10  bearings. 

Pedestal  Data  (1P6E12.5) 

The  program  provides  for  the  option  that  the  pedestals  supporting  the  bearings 
may  be  flexible.  In  that  case,  data  for  the  pedestals  must  be  given  and  NPST 
must  be  set  equal  to  1  (card  2,  Item  3).  If  the  pedestals  are  rigid,  set  NPST 
»  0  and  omit  giving  any  data  for  the  pedestals. 


When  NPST  *  1,  each  bearing  is  supported  in  a  "two-dimensional”  pedestal.  The 
pedestal  is  represented  as  two  separate  masses,  each  mass  on  its  own  spring 
and  dashpot.  The  one  mass-spring-dashpot  system  represents  the  pedestal  char¬ 
acteristics  in  the  x-direction,  (the  vertical  direction)  and  the  other  system 
represents  the  y-direction  (the  horizontal  direction).  There  is  no  coupling 
between  the  two  systems.  In  the  computer  input  there  must  be  one  card  for 
each  rotor  bearing  and  each  card  contains  the  6  items  necessary  to  specify  the 
pedestal  characteristics: 

1  .  The  pedestal  mass  for  the  x-direction,  lbs. 

2 ■  The  pedestal  stiffness  for  the  x-direction,  lbs/inch. 


3  -  The  pedestal  damping  coefficient  for  the  x-direction,  lbs-sec/'inch . 
(Note:  For  the  bearing  films  the  damping  is  given  in  lbs/inch,  whereas  the 
damping  coefficient  in  lbs-sec/inch  is  used  for  the  pedestals). 


5  . 


The  pedestal  mass  for  the  y-direction,  lbs 

The  pedestal  stiffness  for  the  y-direction.  lbs/inch. 
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6  ,  The  pedestal  damping  coefficient  for  the  y-direction,  lbs-sec/inch . 

When  the  bearings  also  nave  angular  stiffness  (i.e>,  NANG  *  1,  item  4,  card  2) 
and  the  pedestals  are  flexible  (i.e.^N^ST  “  1),  the  above  NB  cards  must  be 
followed  by  additional  NB  cards,  one  for  each  bearing,  on  which  are  specified 
the  dynamic  model  of  the  pedestals  for  angular  motion.  Each  card  contains  six 
values : 

1  .  The  pedestal  mass  moment  of  inertia  in  the  x-plane  (around  the  y-axis), 
lbs- inch^ . 

2  ■  The  pedestal  angular  stiffness  in  the  x-plane  (around  the  y-axis),  lbs. 
inch/radian . 


3 .  The  pedestal  angular  damping  coefficient  in  the  x-plane  (around  the  y-axis), 
lbs -inch- sec /radian . 

4  The  pedestal  mass  moment  of  inertia  in  the  y-plane  (around  the  x-axis) 
lbs  - inch^  . 

5  ■  The  pedestal  angular  stiffness  in  the  y-plane  (around  the  x-axis),  lbs. 
inch /radian  . 


b_  The  pedestal  angular  damping  coefficient  in  the  y-plane  (around  the  x-axis), 
lb s- inch- sec / radian  . 


Whirl  Frequency  Ratios  (1P6E12.5) 


When  the  rotor  becomes  unstab Le  it  whirls  in  a  closed  orbit  with  an  angular 
frequency  v  which  normally  is  equal  to  approximately  one  half  the  rotational 
frequency  .0 .  However,  the  exact  value  depends  on  the  rotor  and  the  supporting 
bearings,  and  to  determine  this  the  program  searches  for  the  zero  points  of 
the  stability  determinant.  The  present  input  list  gives  NFR  values  (card  2, 
item  5)  in  sequence  of  the  whirl  frequency  ratio  v/'-O  and  these  values  are  used 
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directly  by  the  program  in  computing  the  stability  determinant.  Thus,  the 
program  is  unable  to  detect  any  possible  solution  outside  the  specified  range. 
Furthermore,  a  zero  point  of  the  determinant  is  only  detected  when  the  determ¬ 
inant  changes  sign  from  one  calculation  to  the  next.  Hence,  if  there  are  two 
tero  points  between  two  consecutive  frequency  ratios,  the  program  is  unable  to 
find  them.  It  is,  therefore,  necessary  to  make  the  increments  sufficiently 
small  (approximately  0.01  or  less).  This  is  most  readily  done  by  giving,  say 
5  or  6  frequency  ratios  in  the  present  input  list,  and  then  specify  an  addi¬ 
tional  subdivision  of  the  intervals  by  means  of  NIT,  item  7,  card  2.  As  an 
example,  the  present  list  may  give  6  values  for  the  frequency  ratio:  0.7, 

0.6,  0.51,  0.49,  0.4,  0.3  (NFR  *6).  In  addition,  NIT  may  be  set  equal  to 
5.  Thereby,  the  stability  determinant  is  computed  at  the  following  26  fre¬ 
quency  ratios:  0.7,  0.68,  0.66,  0.64,  0.62,  0.6,  0.582,  0.564,  0.546,  0.528, 
0.51,  0.506,  0.502,  0.498,  0.494,  0.49,  0.472,  0.454,  0.436,  0.418,  0.40, 

0.38,  0.36,  0.34,  0.32  and  0.30.  Whenever  the  determinant  changes  sign,  its 
value  is  also  computed  at  the  mid-point  of  the  interval  and  these  three  values 
are  then  used  to  compute  the  zero  point  by  quadratic  interpolation.  A  second 
quadratic  interpolation  is  employed  to  get  a  more  accurate  solution. 

Speed  Data  (1P6E12.5) 

This  data  and  the  following  bearing  data  must  be  repeated  NCAL  times  (item  9, 
card  2).  The  speed  data  is  given  on  one  card  with  three  values: 

1  ■  The  initial  speed  of  the  speed  range,  rpm 

2 .  The  final  speed  of  the  speed  range,  rpm 

3 ■  The  speed  increment,  rpm. 

Thus,  the  first  calculation  is  performed  at  a  rotor  speed  equal  to  the  initial 
speed.  Thereafter,  the  speed  is  incremented  by  the  speed  increment,  and  this 
is  repeated  until  reaching  the  final  speed.  The  zero  points  of  the  instability 
determinant  are  found  for  each  speed  and  by  plotting  the  results  as  previously 
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discussed,  the  rotor  speed  can  be  determined  at  which  instability  sets  in.  The 
selected  speed  range  should,  therefore,  be  sufficiently  large  that  it  includes 
the  threshold  speed. 


Eeartng  Data  (8E9.2) 

Each  bearing  is  represented  by  8  dynamic  coefficients  for  radial  motion  such 
that  the  dynamic  bearing  reactions  can  be  expressed  oy  means  of  Eq.  (1)  (or 
more  correctly,  by  Eq.  (5)).  The  8  coefficients  are  given  on  one  card: 


K  K  K  K  u>B  'JjB  'joB  o)B 

xx  xy  yx  yy  xx  xy  yx  yy 


lbs/inch 


where  aj  is  the  angular  speed,  radians/sec,  and  the  8  coefficients  are  computed 
from  the  lubrication  equation  (Reynolds  equation)  as  described  in  volumes  III 
and  VII.  When  the  lubricant  is  incompressible  (INC  =  0,  card  2,  item  6)  the 
coefficients  are  independent  of  whirl  frequency  and  only  one  set  of  values 
should  be  given.  When  the  lubricant  is  compressible  (INC  4  0),  the  coeffici¬ 
ents  are  frequency  dependent  and  one  set  of  coefficients  must  he  given  for 
each  of  the  whirl  frequency  ratios  in  the  input  list,  i.e.  a  total  of  NFR 
cards  (card  2,  item  5).  The  cards  must  be  given  in  the  same  sequence  as  the 
values  in  the  frequency  ratio  list. 


When,  in  the  search  for  the  zero  points  of  the  instability  determinant,  the 
program  performs  calculations  at  frequency  ratios  different  from  the  specified 
values,  the  proper  dynamic  coefficients  are  obtained  by  quadratic  interpolation 
of  the  input  data.  Hence,  a  minimum  of  3  sets  of  coefficients  is  required  when 
INC  4  0,  which  means  NFR  >  3.  When  INC  =  0.  interpolation  is,  of  course,  not 
necessary  and  only  one  set  of  coefficients  is  required. 


If  NANG  =  0  (card  2,  item  4)  the  above  data  is  repeated  for  the  next  hearing 

untiL  the  coefficients  are  specified  for  all  NB  bearings.  If  NANG  4  0,  the 

above  data  is  followed  by  a  similar  set  of  data  of  dynamic  coefficients  for 

angular  motion.  Let  the  dynamic  moment  acting  on  the  rotor  have  the  components 

M  and  M  where  M  is  in  the  x-z  plane  and  M  is  in  the  y-z  plane  (z  is  the 
x  y  x  y 
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coordinate  .long  rh.  rotor  «!.).  The  rotor  a«pUtoda.  .re  x  end  y  end  the  cor- 
responding  rotor  elope,  ere:  6  -  dx/d.  end  »  -  dy/ds.  The  e,o.tlon.  defining 
the  angular  dynamic  coefficients  are: 


dS  d co 

MX  *  “  Gxx  9  '  DXX  dl  '  Gxy  *  *  Dxy  dt 


(34) 


jo  dd> 

My  ■  •  °yx  9  '  V  '  Gyy  ®  "  Bnr  dt 


where  t  1.  tine.  These  equation,  are  completely  .n.logous  to  Eq.  (1).  The 
coefficient,  apply  to  a  given  hearing  geometry,  .  hnoun  static  hearing  o.d  an 

are  function,  of  the  rotor  speed.  For  compressible  lubricants,  the  «.«  fit  » 

tio  such  that  Eq.  (34)  more  properly  should 


also  depend  on  the  whirl  frequency  ra 
be  written: 


M 


Y  e  -  Y  ,  <P 
xx  xy 


(35) 


M 


Y  9  -  Y  <p 
yx  yy 


where : 


Y  =  G  +  i(“)  ui) 

XX  XX  ffl  xx 


(36) 


lmllarly  for  Yjy,  Y^  and  Yyy  Eq.  (35)  1.  analogous  to  Eq.  (5)  and  is  da- 


and  s 

rived  in  the  same  way 


In  the  computer 


input,  the  8  coefficients  are  given  on  one  card: 


G  G  G  G  -sto  Offl  “Dyy  lbs .  inch/radian 

xx  xy  yx  yy  xx  xy  yx  yy 


In  complete  ana 


logy  to  the  input  for  the  radial  dynamic  coefficients,  only  one 
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?ct  of  coefficients  is  required  for  an  incompressible  lubricant  (INC  ■  0)  where¬ 
as  NFK  cards  are  required  for  a  compressible  lubricant  (INC  i*  0). 

When  all  the  required  input,  both  the  radial  coefficients  and  the  angular  coef¬ 
ficients,  has  been  given  for  one  bearing,  the  input  is  repeated  for  the  next 
bearing  until  all  NB  bearings  are  specified. 

Whereas  a  journal  bearing  usually  has  both  radial  and  angular  coefficients,  a 
thrust  bearing  has  only  angular  coefficients  and  its  radial  coefficients  should 
be  set  equal  to  zero.  In  general,  the  angular  coefficients  for  a  Journal  bearing 
are  of  minor  importance  and  may  be  set  equal  to  zero  except  in  unusual  circum¬ 
stances  . 

COMPUTER  OUTPUT 


Referring  to  the  sample  calculation  given  later  it  is  seen  that  the  first  couple 
of  pages  of  computer  output  lists  the  input  data.  Thus,  any  mistake  in  the  input 
data  is  readily  found. 

The  listing  in  the  output  of  the  bearing  data  gives  the  8  dynamic  bearing  coef¬ 
ficients  for  radial  and  angular  motion.  The  angular  coefficients  are  identified 
as  "ANG.KXX",  meaning  G  ,  "ANG.W.BXX",  meaning  a®  ,  and  so  on.  The  first  column 

XX  XX 

in  the  list  gives  the  whirl  frequency  ratio  at  which  the  coefficients  apply.  For 
an  incompressible  lubricant,  the  coefficients  are  independent  of  frequency  but  to 
simplify  the  output  routine  a  value  of  the  frequency  ratio  is  still  given  although 
it  has  no  particular  meaning. 

After  the  listing  of  the  input,  follow  the  results  of  the  calculations.  There  is 
a  list  of  output  for  each  rotor  speed.  The  list  is  preceded  by  a  title  giving 
the  rotor  speed  in  rpm,  and  then  follows  a  four  column  list  of  results.  The  first 
column  gives  the  whirl  frequency  ratio  v/u),  the  second  column  is  the  corresponding 
frequency  ,  radians/sec,  the  third  column  is  the  real  part  of  the  instability  de¬ 
terminant,  and  the  fourth  column  is  the  imaginary  part  of  the  determinant. 


72 


The  purpose  of  the  calculations  is  to  determine  thrie  frequency  ratio  values  at 
which  either  the  real  part  or  the  imaginary  part  of  the  instability  determinant 
are  zero.  The  program  does  this  in  the  following  way:  to  illustrate,  assume 
that  the  real  part  is  positive  at  a  frequency  ratio  of  O.bl  but  becomes  negative 
at  the  next  frequency  ratio  in  the  sequence,  say  0.50.  Then  the  program  cal¬ 
culates  the  determinant  at  the  midpoint  of  the  interval  (i.e.  at  v/co  •  0.505) 
and  through  these  three  values  of  the  determinant  (at  v/co  <*  0.51,  0.505  and  0.50), 
it  passes  a  second  order  polynomial  and  calculates  where  it  becomes  zero.  Let 
this  be  at  v/oj  ■  0.5024230.  This  represents  a  "first  solution".  To  improve  the 
accuracy,  an  additional  calculation  is  performed  with  v/co  **  0.5024230,  and  the 
corresponding  determinant  value  is  used  together  with  two  neighboring  values  to 

compute  a  "second  solution"  which,  then,  is  taken  as  the  final  solution.  Obviously, 

-  6  -8 

the  corresponding  determinant  is  not  exactly  zero  but  it  is  usually  10  to  10 
of  the  neighboring  values.  To  spot  the  solutions  in  the  list  of  output,  simply 
go  through  the  column  of  frequency  ratios  and  each  time  the  results  show  that  an 
Interpolation  takes  place,  the  finally  obtained  frequency  ratio  is  a  solution. 

By  checking  in  the  columns  for  the  real  part  and  the  imaginary  part  of  the  insta¬ 
bility  determinant,  it  is  readily  found  which  part  has  a  zero  point.. 
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INPUT  FORM  FOR  COMPUTER  PROGRAM 


PN4QQ;  THE  THRESHOLD  OF  INSTABILITY  FOR  A  FLEXIBLE  ROTOR  IN  FLUID  FILM  BEARINGS 


Card  1  ( 72  H)  Text 


Card  2  0415) 


1.  NS 

2.  NB 

3 .  NPST  ~ 

n 

4 .  NANG  = 


5.  NFR  - 

6 .  INC  - 

cx 

7 .  NIT 

8 .  METR  = 

9 .  NCAL  » 
10.  INP 


Number  of  rotor  stations  (1  <  NS  <  100) 

Number  of  bearings  (1  <  NB  <  10) 

0:  rigid  pedestals 

1:  flexible  pedestals 

0:  no  angular  bearing  stiffness 

1;  bearings  (and  pedestals  if  NPST  4  0)  have  angular 
stiffness 

Number  of  input  values  of  frequency  ratio  (3  <  NFR  <  50) 

0:  bearing  lubricant  is  incompressible  (liquid) 

1:  bearing  lubricant  is  compressible  (ga3) 

Number  of  subdivisions  of  frequency  ratio  intervals  (1  <  NIT) 
0:  no  diagnostic  METR  ■  0  1:  diagnostic 
Number  of  speed  ranges 

0:  more  input  follows,  starting  with  card  1 
1:  last  set  of  input  data 


Card  3  (1P6E12.5) 


1. 

YM 

Youngs 

modulus 

for  shaft  material,  lbs/inch 

2  ■ 

DENST  = 

Weight 

density 

2 

of  shaft  material,  lbs/inch 

3. 

SHM 

aG,  where  G  is 

2 

shear  modulus,  lbs/inch  ,  and  a  is  shape 

factor  for  shear 
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Rotor  Data  (8E9.2) 


Give  MS  cards  with  7  values  on  each  card: 

1.  Height  at  rotor  station,  lbs. 

2 

2.  Polar  mass  moment  of  inertia  at  rotor  station,  lbs. inch 

2 

3.  Transverse  mass  moment  of  inertia  at  rotor  station,  lbs. inch 

4.  Length  of  shaft  section  to  next  station,  inch 

5.  Outer  shaft  diameter  for  cross-sectional  moment  of  inertia,  inch 

6.  Outer  shaft  diameter  for  shaft  mass,  inch 

7.  Inner  shaft  diameter,  inch 


Bearing  Stations  (1415) 

List  the  NB  rotor  stations  at  which  there  are  bearings 


Pedestal  Data  (1P6E12.5) 

This  data  only  applies  when  NPST  +  0  (card  2,  item  3).  Give  NB  cards  with  6 
values  per  card,  one  card  per  bearing: 

1.  Pedestal  vibratory  mass,  x-direction,  lbs. 

2.  Pedestal  stiffness,  x-direction,  lbs/inch 

3.  Pedestal  damping  coefficient,  x-direction,  lbs. sec/inch 

4.  Pedestal  vibratory  mass,  y-direction,  lbs. 

5.  Pedestal  stiffness,  y-direction,  lbs/inch 

6.  Pedestal  damping  coefficient,  y-direction,  lbs. sec/inch 


If  also  NANG  /  0  (card  2,  item  4),  give  additional  NB  cards,  6  values  per  card 

2 

1.  Pedestal  vibratory  mass  moment  of  inertia,  x-plane,  lbs. inch 

2.  Pedestal  angular  stiffness,  x-plane,  lbs .  inch/radian 

3.  Pedestal  angular  damping  coefficient,  x-plane,  lbs . inch . sec/radian 

2 

4.  Fedestal  vibratory  mass  moment  of  iaartia,  y-plane,  lbs. inch 

5.  Pedestal  angular  stiffness,  y-plane,  1 bs  .  inch/radian 

6.  Pedestal  angular  damping  coefficient,  y-plane,  lbs .inch. sec/radian 
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Whirl  Frequency  Ratios  (1P6E12..5}. 


List  NFR  values  of  the  whirl  frequency 
Note:  The  remaining  data  must  be  repeated  NCAL  times 

Speed  Data  (1P6E12.5) 

1.  Initial  speed,  rpm 

2.  Final  speed,  rpm 

3.  Speed  increment,  rpm 


ratio  in  sequence,  6  values  per  card 


Be  a  rim:  Data  (8E9.21 


,  i  it  rwr  =  Q  with  8  dynamic  translatory 
ji  0,  or  l  card  only  it  INC  =  0,  wren  / 


Give  NFR  cards  if  INC 
coefficients  per  card: 

Kxy  v  Sy  "*»  “V  lbs/t”ch 

I£  also  ,  o,  *U  1.  *  «  «  *  —  8  d>"“1C  a"8U'ar 

cottf >c rents  per  card: 


r  a  G  ojD  cuD  'JoD 
xx  xy  > 


G  L’-x  uyy  “xx  -xy  “yx  "  yy 


o£)  lbs.  inch/radian 


Thcie  is  one  se 
repeated  NB  times 


,  ^  p  the  above  bearing  data  is 

t  of  coefficients  per  bearing,  i.e.  tn 
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DIMENSION  mm  (  nO)t**IP(  80 )  •  H 1 1  (  rtO)  .KL(  .80)  *HS(  80l»«*(  80>» 

C  PN400*STAftlLITf  ImkEshOLU  OF  FUfcAlbLfc  KOTDK 

1RD(  60>«MJ(  btU.OvM  rtO )  (DMA  (  oO).UM8<  SO)  «Al  (  B0>*A2(  80)  * 

2A3(  e0)*A8(  Huli^l  111)  «AM  bO)*AM  aO>*AM(  60)*BN(  8G)»l)N(  80>» 
3FH01  (80)  .XL (4)  .SL  (8)  *«L  (**)  ♦  JL  (8)  *XK(8>  *SM(8)  tab (81  *VH(8>  *LB(  10) 
DIMENSION  MMX  (2»  HD  «mi)A  (2*  10)  »PMY  (2*  10)  *P*Y  (2*10)  • 

1P0Y(2«10>  *sa<?*») 

DIMENSION  uC(m*80*10)  «i)A(b. 80.10)  *  o  V  (2.8.10) 

COMMON  AC  (4.8)  .AS (8.8)  ♦1)Um)TS 

200  WEAIK5*  1  «»0 ) 

Hf  AO  (5*101)  NS.N*-.  .hsi  .r-Ai'  bi«fK«  Inc  * N I T * ME T« *NLAL»  IM* 

PF  AD (5. 1 0? )  Y*  .ntNbl  •'•hm 
*)MITt(tt«10()) 

AMITE (b*lui) 

#MITt'(b.  l'jHjNStlMtl.I.Pa  I  .NAMi.NFN*  INC f  N IT  ♦  MET M* NCAL*  IMP 
KH(Tt(8. 10b) 

4W  i  Tt  (6 « 1  On )  Y  *  ••it  Nat  *  jn*> 

WKITF  (6.107) 

VfPITt  (o«  I'M) 

DO  201  J*l.Nb 

MFALKS*  108)  r.n  (  J)  IM  (  J)  «wl  I  ( 0)  ♦ML  ( J)  *KS(  J)  * N*  ( J)  »MU  ( J ) 

201  HHJ1F  <6. 1 10)  J.-'-U  J)  ...  )M  ( j)  «H  l  |  ( j)  ,«(_( J)  »MS(  J)  «Kl»(Jl  *HO(  J) 

PF  Al<  (5*  1 0 1 )  (I  ( J)  •  0=  I  »N*) 

MMtTF.  (6. 1 1 1 ) 

WRITE  (6*  10  i  >  (Lr*(o)  •  J=l  «.*n) 

202  AM<j=3Sb.068 

AMSL=0. 0001)01 

208  MS  (.'IS)  =1.0 
DF.NST  =i'r.NST/'  Af’> 

00  210  0=1. .ib 
pM(J)=HM(  J)  /A.'Yt 
H1P( J)  =X|U(J)  /A  <5 
RIT(J)=8lT(  Jl/ws 
C2=mik  j) 

C3=KS( J) 

HS  ( J)  =0.  uv/ilh  7  =  «  t  "»  (  (LJ««8)  - (C2**8)  ) 

C4  =  K«  <  J) 

«W(  J)  =0.  /rts.)8-4n»  <1. <♦<»!. s -!.<:*(.<;)  “OtMSl 
IF(C8)  20a.20S.2on 

205  PJ(o)=u.U 
GO  TO  2 f»7 

206  HJ(J)  =  (C8*C8*CV*t£>/ 12.0 

207  C?=1.5  7o78tj*5n«MC3’>('J-C2ftC2) 

I F  ( L2 )  206.20 7.202 

204  PD( J) =  nS( J) /C2 
GO  TO  2 1 o 

209  Ml  i  ( J )  =0  •  u 

210  CONTI  wot 

IF(NMSI)  211. 215.211 

211  WRIT!  (6.112) 

WSITfc (6.1 13) 

DO  212  0*1.. *8 

RE Au (b« 102)  Pt-  a  ( 1  ♦  J)  .HN A  (  l  •  j)  .MDA  (  1  .  J)  *MMY  (1.  J)  *PKY  ( l.  J)  ♦ 

1  MO Y ( 1 . J ) 

WWITt  (ft*  1  10)l.n<  J)  *MMA  (  1  «  J)  .  PK  A  (  l.'j)  *  MDX  ( 1  »  J)  .PMY(1»0>  .PNY  (1*  J)  . 
1MDY (1.0) 

PMX  (  1  .  J)  sui«x  (  J  ,  J)  /AMJj 

212  PPY(l.0)=PMY(l*0)/»Mb 
IF  (NANG)  213.215.213 

213  wMITF (6.1 1 8 ) 

DO  218  0=1.  imm 

HE AD (5. 102)  PMX (2 .0) .PK A (2,0 )  .POA  (2.  J)  * PMY  (2.0 >  * PKY  (2*0)  * 

lPny (2.0) 

WkITtr  (b.  1  1  0  )  Lu  (  0)  .MMX  (2 .0)  ,MKX  (2.0)  * PDX  (2.0)  .PMY(2.0>  .PKY  (2.0)  . 

1 POY (2.0) 

P’-iA  (2  •  0)  =mmX  (2.  J)  /A».b 
218  PMY (2*0) =PMY (2. J) 
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215  REAO<5.102)  (FhOI Jj) J J=l, \FR) 

WR!TE(6.116> 

WRITE(o.lOft)  (KkUI  (j)  «J=1.|\Fh) 

IC*1  70 

301  REALMS.  1U2)  SPsT  .5rF  I'i.ShnC 

K1«NFW  72 

IF (INC)  303. 302.301  73 

302  K 1= 1  74 

303  xWIlK (6.11ft) 

WRlTtlb.lOblSpsT .bPrw.bKNC 
wRlTtl4.ll 7) 

00  310  KM..14  78 

WR!TE(6.1l6)L4<«) 

w«ITt(6.11*) 

00  J06  J=l.Kl  til 

RFA0(5.109)  (DC(I.J.iU.l  =  l,8> 

306  *R?  f  t  (  6 ♦  1 2  l  )  F w 1  {  J)  ,  (*jC(  l.J.M  .1*1.6) 

IF(NANO)  3u 7*310* 30 7  84 

307  ..UITt  (6*120) 

00  J0«  J=1.M  86 

HFAlMS.  109)  (•>•>(  I  .  l.M  .  l=l.d) 

309  WK1TF  (o.  121  )F-.Jl  (  j ).(.;«( 1 ,  j.ft)  .  1  =  )  ,h) 

310  CONTlNur.  89 

SPI)  =  5PSI  90 

320  ANSP=0. 104  71<.7b»b-*t;  91 

WRITt (6.122) 5Rn 
WR  II  t  ( 6 . 1 2  3 ) 

10=1  94 

FRm  =  KR*>1 ( 1 )  95 

KR=0  96 

KE=0  97 

K 1 =1  98 

K?  =  0  99 

K3=0  100 

0F-*sl9[T  101 

DFR* (FwO 1  (2)  -F  •.»)  /i>  *  102 

325  FR0=Am5P*Fh*  103 

FR02=FHvJ*Fw«j  104 

00  3/6  J*I...5  105 

0VA( J)=FwO/*RV(j)  106 

O^A  ( J)  =-Fr(  j2*vl  T  (..)  107 

326  DMR(J)=FM)*ANSP*Hlr I J)  108 

330  00  375  J=l.'iR  109 

IF(INC)  3)2. 331.3)2  110 

331  11=1  111 

GO  TO  3  112 

332  IF(r?)  lift.  j.3J.  .j-<  113 

333  Il=Iu  114 

334  00  3 ji'  1*1.6  115 

SM(l.I)*OC(t.n.J)  116 

IF  (fJiiNO)  iJj.Jjft.  117 

335  $r.  (2.  I  )  =CA(  I  ,  i  i  .  j)  118 

33ft  C0.9Tli.ur.  119 

GO  In  120 

339  1 1  *  1 u  121 

I  ?  =  2  122 

1.3=0  123 

IF  ( ft 3 )  340.J-4li.  jj-.  124 

339  11*19*1  125 

340  IF  (..FP-  I  1  )  ISO.  15b.  341  126 

341  C4  =  Fk01  (ID  127 

Cl=FWul ( 11*1 )*C4  128 

C2=C4-Fh>u1  ( 1 1  - 1  )  129 

C J=C 1*02  130 

C9=Fkw-C4  131 

DO  355  1-1.8  132 

C5*0C(I.I1.J)  133 

C6= ( DC ( I • 1 1 ♦ I . J) -C5) / (C1*C3)  134 

C4=(0C( 1. 1 1-1 . J)-Ca)/ lC2*CJ)  135 
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C7«C6*C8 

£6*C2*C6-Cl*C8 

C8*CS*C9*(C6*C9*C7) 

IF(I1-I2)  343.342*343 

342  se<i*n*c8 

IF (NANG)  347*355.347 

343  IF (13)  344*344 • 345 

344  SB  ( I  •  I )  *C8/C2*  (FH#-FK'Jl  ( 1 1-1 ) ) 

IF(NANG)  ^47.355*347 

345  SB(1.I)=SB(1*I)*C8/C1*(FKG1  (I1*1)-FF<«> 

IF (NANG)  347*355*347 

347  C5=QA ( I « 1 1 « J) 

C6*(DA(l,ll*l.J)-Cb)/(Cl*C3) 

C9* (UA (1*11-1. J) “05) / (02*03) 

C7=C6*Crt 

C6*C2*C6-C1*C8 

C8=C5*C9*(C6*C9*C/) 

IF(I1-I2)  349. 346.34* 

348  SB (2. I) =08 
GO  TO  355 

349  IF  (13)  350.350*352 

350  SB(2.I)=C8/C2*(FK<i-FKJU11-1)  ) 

GO  TO  355 

352  SB(2.I)=SB(2«I)  *08/01*  (FKt*l  (11*1  )-F«F) 

355  CONTINUE 

356  IF  ( *'  3)  359.357. 354 

357  IF ( I 1-12)  359.359*358 

358  11=11-1 
I2=NFk-1 
13=1 

GO  TO  34 1 

359  DO  3b0  i=l«H 
0V(1*I.J)=0.0 

360  DV(2*I*J)=0.0 
11  =  1 

361  oo  362  1=5.8 

362  SB  ( 1 1 «  1 )  *FH  »*Sr*  ( 1 1  ♦  1 ) 

IF(NPST)  365*363.366 

363  00  364  1=1*8 

364  0V(Il.I.J)=5B(n.l) 

GO  TO  370 

365  C1=MK*  ( 1 1*  J)-F.«1c*(-4a  ( 11  *J) 

C2=FNO*8GAdl.J) 

C1=PKY( I l.JJ-F 802*84 Y( tl.j) 

C4=Fft(J*r’0Y  <  II  .  J) 

C5=Sd»l 1*1) »C1 
C6=Sh( 1 1 .5) *C2 
C7  =  SM  I 1 .4) *C  3 
C8  =  SH ( 1 1 «  B ) *C4 

01 =05*07-06*06 -Sh  ( !  1  *2)  *5b  ( II  *3)  *S8(  i  1  *6)  *SB(  1 1  *7) 
D2=C5*C'<*C6*C7-Sm ( I  I , 2) *58 ( II  *  7 ) -SH ( 1 1  *  3) *SH (11*6) 
IF  (At»S(Ol  >-AHS<02>  I  Jr>7  *  366  *  36b 

366  03=02/01 

04  =  01*0  )*02 
01=1.0/04 
02=-03*0l 
GO  TO  368 
347  03=01/02 

04=02*03*01 

D?=-l .0/04 

01 =-03*02 

368  03=C1*OI-C2*02 
D4=C1*D2*C2*01 
05=1 .0-03tC/»D4*C8 
06=-03*CH-U4*C7 
07  (I  l»l*J)=C’.  *05-02*06 
OV  (11»5»J)=C1  *IJfc*C2*05 
05=C3*5H  (11*2)  -C4*'3ri  (11*6) 

06«C4*S8 (11,2) *C3*S8 (11,6) 
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DV ( 1 1 *2* j)  =f)3*05-GH*i>o  205 

DVlI1.6.j)=03*06*U4*U5  206 

DSxC3*SM<11.3)-C4*->h<11.7)  207 

06=C4«:>-t<  U  .31  »C 3*3"  it  1  •  / 1  208 

DV(n»3*j)=03*u5-')‘*«j-->  209 

0VU1.7.jI=o3*06*'J4»u:>  210 

D3=C3*01-C4«t)2  211 

04=C3*02*C4*Ul  212 

05=1 .0-DJ*C-3«n<4*C6  213 

06=-03*C6-04*C:>  21* 

DV<ll.4.j)=CJ*05-C440o  215 

OV  (ll*B»J)=C3®L)6*C1**l)'j  216 

370  IF(NANG)  371* 3 75*371  217 

371  IF(Il-l)  372*372* i/5  218 

372  11=2  219 

GO  10  JM  220 

375  CONTlNtJt  221 

IF(MtlH)  3rtl. 352,3*2  222 

3fl  I1=Iams(-F1h)  22 3 

rfWIIr  lo.l2l IF**.  (>v< 1 1 .1.1) . 1=1.3 

GO  Tu  4»o  225 

332  COf.'T  I  Nut  226 

1 1  =Ni“ 1  227 

00  420  J= 1 . 1 1  228 

C1=MS<J)  229 

C2=FH02*nk( J)  230 

C3=-7(XJ)  231 

C10=*JU)  232 

C4=C2/C1  233 

CS=S<JHf(t4)  23* 

C6=SuwT(C5)  235 

C7=f<L<J>  236 

IF  <Cb*C  7-u.u J)  4n.4ii.4U  237 

411  C3=C2*C7  238 

Dl=C7*C7*l>  239 

02=C7*l)  1  2*0 

D3=C  7*02  2*1 

D4  =  r)l*<C3-ClO)  2*2 

D5=C  7/C 1  2*3 

07’  =  C7«Uo/2.0  2** 

D7  =  C  7®Uo/ 3 .  j  2*5 

DA  =  C3-*C  1 0  2*6 

0F=tJH*l)3*CJ*Cli)  2*7 

A  1  I  J)  =1  .0*1)  3/24.0  248 

A2  (  J)  =<>  l  <  J)  “!'■*  2*9 

A3  ( J)  =  t  1 .0  —  >4/  i.o*OJ/  l/U.onc/  250 

A4  (  J )  =IV/4. 0  251 

AS(j)=C1«A4< j)  252 

AM  J)  =C2«A  J  ( J)  253 

lMJ|s(1.0*7,'i*C4»U*(  10H.J/12U. 01*05  25* 

4‘  1  J)  =l.u-')4/“>.n*-i3/3f*U.O  255 

A7  ;  J)  =  C->*~in  <  j) /£-.  f-*c.  /  256 

riN  (  J  )  =  llo*M<-  <  J  )  257 

l)f:  (  J)  =  <  1 . 0*4.n'*r«*4i/M»U  1-  <CJ-ClU)  *2.0*05  258 

GO  Tc  -.✓u  259 

412  r,-=t-."  i  u  i*ri .  )  e»2i  260 

c^=rs*  ttj-ciM  261 

1  F  j  .  u  III!  2 )  1  1  J  •  4  1  j  •  H  1  •*  262 

4l.j  Cl  1  =  1  .u*k. /-**'.•'  263 

C  J  =  1  ,  0  *  u  .  264 

GO  lo  -*15  265 

414  C*=SvkT  <  1  *u*C->  266 

Cl  1  (Crt)  267 

415  05=C5»  <  C-i-l  i 1  268 

00=05“ <C-*<  X  2t9 

.  D9  =  05*t)o  270 

IF  (6w5tC**/Cf1)”'J.00  02>  410*416*417  271 
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416  C1»*6-S/C«*C9  272 

O3aCr>Cll*a.»i-0tt»  273 

04*C6*C11*(1.U*PB)  274 

GO  TO  4 1 «?  276 

417  DSsSC’ftKCS)  276 

04*SU«T<06)  277 

413  07*03*05  278 

08*04*04  279 

01*03*07  280 

02*04*07  281 

VC*COS( 02 1/09  282 

VS=SlN<D2>/09  283 

8MC*f*O(01)  284 

BHS*l.«/fOC  285 

Ol*(».5*leKC*6MSJ/,J,»  286 

02*0.5*  (o*"C-HMb )  /l“»  287 

A1 O) *04*01 *05«VL  288 

A2  (J)  *05*01  *li4»v(.  289 

A3(  J)*03*02*U4*V5  290 

A6 ( J) =C2*AJ< J)  291 

BN< J)=(01-vC)/Cl  292 

A7(J)*C2*(Ul-wCl  293 

A4(J)*C5*<U4*n2-L.l*vS)  294 

A5(J)*C1*A4(J)  295 

Cfl*Cl*C5  296 

AM J)  =  (Oe*d2*i>7*vS)/(  >  297 

0N(  J}*CJ7*02-:>r1*V5)/C<-  298 

420  CONTINUf  299 

00  470  1=1.4  300 

00  4b2  J= 1.4  301 

XL(J)=O.0  302 

SL(J)=0.0  303 

BL  ( J) =0.0  304 

452  VL(J)=0.0  305 

11*1*1-1  306 

1 F  ( 1-2 )  4‘aJ, 4-3.1. 454  307 

453  XL«I1)=»4SL  308 

GO  TO  455  309 

454  11=11-4  310 

SL  ( 1 1 )  *A*«SL  311 

455  11*1  312 

I2=LM1)  313 

00  465  J*1.N5  314 

00  456  L=1 .4  315 

HP(L)=Hl_(L)*04a(j)»t.L(L)  316 

VR(L)*VL(U  *‘)VA(j)*/LlL)  317 

XR(L)=XL(L)  318 

456  SR(L)=SL(L)  319 

BR(l)=HR(l)-0‘(tMj>«SL(4>  320 

a»(?)=«'-((2)  *0'*i  ( j )  *Sl.  (  J )  321 

BR(3)  *e«(  31  *0MrJ  ( J)  *SL  (2)  322 

3R(4)=rrM4)-D64(j)*SL.(n  323 

IF ( 1 2- J)  46  0 .457 .4o0  324 

457  8R(1)=b*(1)*SL(1)*0V(<'.I*I1)-5L(2)*i)V(2.S»I1J  *SL  (3)  *0912*2. 1 1 )  —  325 

1SL  (4)  *0V  (2.6.  ID  326 

dR(2)*HK(2)*5L<l)*OV(2.S,ll)*3L(2)*UV<2,l.Il)*St(3)*09<2*6»H)*  327 

1SL (4) *00 (2.2.1 1 )  328 

BP ( 3) =8R ( 3 > *SL (1) *OV (2.3. 1 l)-5L( 2) *0V (2.7. 1 1)*SL  (31 *0V12*4» III-  329 

lSL(4)*0V(2.«.in  330 

BR  (4 1  =RR  (4 )  *SL  ( 1 1  *i'V  ( 2. 7 . 1 1 1  ♦  5L  (2)  *OV  (2.3. 1 1 1  *SL  ( 3) *0V(2 .8.  Ill*  331 

1SL(4)*0v(2.4,I1>  332 

VW(1>*VR(1)-XL(1)*UV(1»1»11)*al(2)  *OV  (1.5. Ill -XL  ( 3) *OV(l»2» I 11*  333 

1XL ( 4 ) *0v (1.6.11)  334 

VR ( 2 1  * VK ( 2) -XL ( 1 1  *OV ( 1 ,5  » 1 1 1 -AL ( 2) *0V ( 1 » 1 . 1 1 ) -AL ( 3) *0¥tl  »6»  1 11-  335 

1XL(4)*0V(1.2*I1)  336 

VR(3)=VR(3)-XL(l)*0V(l»3.Il)*XL(2)*0V(l«7»ll)-XL(3) *0¥(1»4» 1 1)  ♦  337 

1  XL ( 4 1 *U v ( 1 ♦ d. 1 1 1  338 

VR (4) =VR (4) -XL ( 1 1 *OV ( 1 . 7  *  1 1 1 -XL (2) *0V  1 1 »  3. 1 1 1 -XL (3)*0K11*8.I11—  339 

lXL(4)*DV(1.4,Ii)  340 
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n«n*i  34i 

IF(Nh-II)  458*464 ***59  342 

458  |;?»NW  343 

60  To  460  344 

454  I  ?»Lk'  fill  345 

460  IFTNS-J)  466*666. •♦ni  346 

46 1  00  46?  L  *  l  ♦•»  347 

XL  <L)  «A?  ( J)  «  x~  (L  )  ♦  ;»  J(  J)  <*Sh  <L>  *rN(  J)*8W<L)  ♦DMJ)«V«(L>  348 

SL  <L»«A4  t  J1  <l.)  *«1  (  J)  *»b6  (L)  *A*t  J)  ♦8N1J)*VH(L>  349 

BL  <L>  »A7  ( J)  *x~  <L>  ♦A-'  (  i)  <*68  <L)  *Al  t  J)*BHML>  *A3  ( J)  (L)  350 

46?  VL  <L)«A6(.M  *«- <L)  ♦<*  /IJ)  #iHU  *m4<J)*88(L)  *A21 J) *VH It)  35i 

465  CONTINUE  352 

AC(l.I)sv-'U)  353 

AC(?»l)*v<*(  j)  354 

AC  (3* I ) *6« ( 1 )  355 

AC  (4»  I )  =*)H  (  s)  356 

As  ( 1 .  1 1  =*•*  (?)  357 

AS (?• I ) =v~ (4)  358 

AS ( 3. J )  =ri'  I  ?  1  359 

AS (4, 1 ) =64 (4)  360 

470  CONT | 361 
CALL  CUtl* 

*iW|TF  (  6  •  1 IJ6  )  f  •  0  ( t  •  O  1  3 

364 

480  IF(M-)  6ii  »bt->  365 

4(5 1  IF  (K  1  )  4(  J,/*MJ.4S?  366 

<*»2  Rl  =  0  367 

DC l=b 1C  368 

TCl=n-*  369 

1*0  T(;  -H6  370 

46 3  IF  (0tC«  »C  1  )  446  M?t  •-*-•*  371 

484  l)C  1  =I.)T C  372 

TCl=FH*.  373 

486  I  F  ( f/TS“OS  1  )  6i/3  *460  •4.-»b  374 

466  L)Sl='lTb  375 

TSI=(h»  376 

IF(NK*-IU)  jtO  •  65<>  *4n  /  377 

487  IF  (Nl  T-l -X  J)  4rt4<46<»«480  378 

469  K3=63*l  379 

K<?  =  1  380 

FFW=M*»*QFX  381 

60  To  j?3  382 

4M4  K  i=ii  383 

K?=o  384 

IO=lo*l  385 

F ;■>►=(•  HU  <  Io>  386 

IF  (•<*•  w- (  (,  >  4Hf«4V0»44J  387 

600  0Fw=G.l  388 

60  Tii  j?6  389 

441  i)F*,  =  r«Il  390 

l'F*=  (F  h  jl  ( 1  j*  1 )  -lh  4 > /Lih  n  391 

60  TO  ,5f6  392 

44?  Kwrl  393 

•.)C’5=|/1  C  394 

If  <  =  FF  *  395 

OS  1SHI6  3*>6 

rsi=fH«  397 

FH-=(fH«*TCl)/2.U  398 

x?=l  399 

C-i  i  lu  J?S  400 

445  KC=*1  401 

X  1  =  TC 1  402 

X  ?  =  F  HW  403 

x  3  =  TC3  404 

r 1 =00 1  405 

Y?  =  l>Tc  406 

Y3=t»Cj  407 

GO  TO  5?0  408 

446  IF<KH*1)  601 .447*  467  409 
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♦97  KR«-2 

♦  10 

IF<FRlr-X21  499#498»49a 

♦  11 

♦96  X 1 »X2 

♦  12 

Y1*Y2 

♦  13 

GO  TO  500 

414 

♦99  X3*X2 

♦  15 

Y3»Y2 

♦  16 

500  X2=FHW 

♦  17 

Y2*0TC 

♦  18 

GO  TO  520 

♦  19 

501  KR*0 

♦20 

DTS=0S3 

♦21 

F9W= TS3 

♦22 

DC  1 =DC3 

♦23 

TC1=TC3 

424 

GO  TO  4*5 

♦25 

505  KF  =  1 

♦26 

0S3=0TS 

♦27 

TS3=FWw 

♦28 

FHW=(Fkx*TS1)/2.0 

♦29 

K?=l 

♦30 

GO  TO  325 

♦31 

506  KE=~1 

♦32 

Xl=TSl 

♦33 

X2=F«W 

♦34 

X3=TS3 

♦35 

Y1=0S1 

♦36 

Y2=0TS 

♦37 

Y3=DS3 

♦38 

GO  TO  520 

♦39 

507  IFIKE+l)  512»508*50d 

♦♦0 

508  K(-=-2 

♦♦1 

IF(Fww-X2)  51  0  «50'*»60'» 

♦♦2 

509  X  1  =  X2 

♦♦3 

Y1=Y2 

♦♦♦ 

GO  TO  511 

♦♦5 

510  X3-X2 

♦♦6 

Y3  =  Y2 

♦♦7 

511  XC=FWw 

♦♦8 

Y2=0TS 

♦♦9 

GO  TO  520 

♦50 

512  KF  =  0 

♦51 

DT  S=i)S3 

♦52 

f  Pw  =  T S3 

♦53 

GO  TO  4«* 

♦54 

520  CI=*J-'2 

455 

C2  =  *  2-x 1 

♦56 

C 3  =  x  3-x  1 

♦57 

C4  =  (  Y3-Y2)  /  (Cl*C3) 

♦  58 

0=<Yl-Ye)/<l.2*Cj) 

♦59 

C5  =  C2*CH-Cl'M;n 

♦60 

C4=C4  *  C6 

♦  61 

I F  (  C**  )  5c2»521*522 

♦62 

521  C6=-Y?/C5 

♦63 

GO  10  525 

♦64 

522  05  =  0.5/04*0 

♦  6b 

C4  =  4nS(C5*C:j-r2/C4) 

466 

C4=50kT (C«) 

467 

IF  (C5)  523  *524 « 524 

468 

523  C4=-C4 

469 

524  C5=0-C5 

470 

525  FG'v-/2*C^ 

471 

GO  TO  32o 

472 

550  S P 0  =  S K 0 ♦ 5 2 1 J C 

473 

IF  (SPHm  +  O.  00001  -akiJ)  t>52  » 552  »  32 0 

474 

5S2  IF(nCAL-IC)  5-d4,554»5->3 

475 

553  I C= I C  *  1 

476 

GO  ’0  J  )  1 

477 

554  IF(tNP)  t5d»2U0*5'35 

478 
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5S5  CALL  FAIT 

100  FORMAT (/ eh  480 

1  )  481 

101  FORMAT(l<,l3)  482 

102  F0R'-iAT(bhl2.b) 

103  FORMAT  MnHiIST  A  f  IGN;>  .jO.nFub  f-Eu.FLEA  ANG.FLFX  .MO.FRLU  INCOMPK  484 

1  SIMP  I  VS  M*  T-MC  SPttOS  I*Rul>  488 

104  FORMAT  ( Ib.RlR)  486 

105  FOKMATUjHO  Y'.>uNL'b  M  Ji  -  •  5  AfSnuE  1*5  1  TY  (  bMfc  Arr  FALT>*G>  48? 

106  FORMAT  (Ht  1  *♦ . *3 J 

107  FORMA  T  (  I  1mu~uT  .)*  i.«Ia)  489 

106  FORMAT  UO*mo->T A!  Iui,  r,«SS  POLAR  MOM.  IN.  iHANbV.MOll.  IN  L  490 

IFNOTn  OT.i'UI'j!  IF  h  )  will  .01 A  (MASS)  lN.MtR  oiA.)  491 

109  FORMAT Mtv. 2)  492 

110  FORMA  T  (  |b  *F  1  N.r,  ,  'ft-  1  A  .rv) 

111  format  1 1  VHuRf- ARi.\f*  stations)  494 

112  For1-  uT  (  1  -omr'F  1  S  1  «l  !/-l  A)  495 

113  H>r.'-»TI  r<)hu-«T-1  1 1* •.  -»SS-A  STIFFNESS"*  DAMPING"*  MA  496 

1SS-Y  OI^V  >y*»  l  AMR  li.G"  Y )  497 

114  FORMAT  t  yii*-0->T  M  lOS  r-i-n.  INtRT-A  ANG.STIFFN-*  ANG.DAMF'-*  MOM.  498 

1INFkT-Y  **>»•>!  IFF n-t  aNw.oanr-t)  499 

115  FORMAT (2jmuTxJ AL  F  KfcojtNLY  hATIO^)  500 

lib  FORMA  I  (•*£*•  i  Ill  ft*L  Sri- tO  M  l«L  SPt.t.0  SPEEO  INCR.)  501 

117  FORMAT  (/R'f-l'tfwi.-'IC  OrMRlNc  CutFr ICIENTb)  502 

116  FORMAT  (  /cirij-ii  !•-'  i  -•j  At  ro  I  jn  -jUU  JM  J)  503 

119  FORMAT  (  1  i-  11  i. A  I  i  cnAjrKxxi  0*  JrtKA  T10XJMKY*  1  OX  3FIKY  Y9X5HM*bXX8X5  504 

IHV.  *««Y»<:A.'i  0<"»->n(ii'.-*Y)  505 

120  FOrRmT  (  ljHil  F  -F  i,  .  Ru  I  l  Or  a  ihiiN'j  .  N  A AoA  7HA..G  . K  A YbA  7HANG.KYXbX7HANG.KY  506 

1  Y5»  ANb  .  wteri  Y  X4XRHANG • Y Y )  507 

121  F  ORMR  I  I2fcl4.5*/*1  .'•“>) 

122  FoRFi aT  ( ///  1  3-.- «  (1  .. r<  irttU^tUiOinH  RRm ) 

123  FORMAT  (SbMil  F«Fij.«a;  tu  FKtuutNCY  Rt(GETtRM)  IMDETERM))  510 

End  511 

SUrlwIJUIlNt  CPF.Ir 

C  SllRROUllNr  F  JR  RN400 

OI.ifnSIun  Imi,<(a.j) 

CnvwOFj  AC  (<••<*)  «AS(**«4>  .UlL.Olb 

C  INITIALIZATION  4 

N=4 

OTC=l  .o 
OTS  =  ti.ll 

00  10  J=1»N  6 

10  T NO A ( J*  J ) =  J  7 

C  CRFCK  FOR  2£r<>  mATkIa  8 

OO  15  1=1. N  9 

I«*=l  10 

IC.L  =  1  11 

DO  14  J=1.N  12 

IF  (AHS  («C  <  I  .  J)  )  *A--»(A>(  1  .  j»>)  12*12*11  13 

11  I Rw  =  0  14 

12  IF"  (APS  (AC  (J*  1  )  )  *Ar>b<A->(J»  i)  )  I  14*14*13  15 

13  I  CL  =v  16 

14  CON'T  I NUF  17 

IF(|RR*ILL>  l“>*lb*4b  18 

1 5  CONT  1  NOr.  19 

C  SF AkCm  FGR  PlvoT  tLFMtNT  20 

00  4  i  |  =  1.  A:  21 

Ts  T  =  b . 0  22 

Oo  24  .1=1  •  n  23 

IFIIM'*(j.jl-l)  1/.24.I7  24 

17  00  «  =  l-r:  25 

IF ( INOX <*.3>-lJ  lb.23.4b  26 

1H  Cl=SuRT(4C<J.K>*AClJ.K>*AS<J*f\>*AS(J*Kn  27 

IF I TST-C1 )  22.23.23  28 

22  I AR= j  29 

I Cl -x  30 

TS  T  =C  1  31 

23  CONTINUE  32 
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24  continue  33 

IN0X<ICL*3)*INUX(ICL*3)*1  34 

lNOX(I*l)»IRN  35 

INOX  1 1 *2 ) * ICL  36 

C  INTERCHANGE  RO*S  lu  RUl  HlvOT  ELEMENT  ON  OIA6UNAL  37 

IFIIRX-ICL)  25*29*25 

25  OTC»-OTC 
DTS«-OTS 

DO  26  L* I *N 

CI*AC  ( IRw.L)  40 

C2=AS<  H».L)  41 

AC(IRw*L)=AC(ICL*L)  42 

AS(IRW*L)=AS<ICL*L)  43 

AC ( ICL«L) =C1  44 

26  AS ( 1CL*L ) 3C2  45 

C  DIVIDE  RIvof  RO»  yr  PIVOT  ELEMtM  46 

29  C1=AC ( ICL* ICL )  47 

C2*AS ( ICL* ICL)  48 

TST*OTC 

0TC*C1*TST-C2*DTS 

0TS*C2*TST*C1»JTS 

AC ( ICL* ICL) =1 .0  49 

AS  ( ICL  *  ICL )  -0  <0  50 

IF ( *BS (Cl ) ~A4S (C2) )  31.30.30  51 

30  TST=C2/C1  52 

C1=1.0/(C1*TST*C2)  53 

C?=TST*Cl  54 

GO  TO  32  55 

31  TST=C1/C2  56 

C?=l .0/ (C2*TST*C1 )  57 

C1=TST  *C2  58 

32  DO  33  L=1.N  59 

TST=AC< 1CL.L)  60 

AC( ICL*L)=Cl*TST*C2*Ab< ICL.LJ  61 

33  AS(ICL.L)=C1*AS(ICL.L)-C2*I5T  62 

C  REDUCE  NON-PIVOT  ROWS  63 

36  00  <*1  L 1 3 1  •  N  64 

IF (Ll-ICL)  37,41,3/  65 

37  ClsAClLl.ICL)  66 

C2=AS ( L 1  * ICL)  67 

AC(L1*ICL)30.0  68 

ASlLl.ICLJsO.O  69 

DO  3*  L  =  1»N  70 

AC(L1.L)=AC(L1.L)-C1*aC(KL,L)*C2*AS< ICL,L)  71 

38  AS(C1*L)=AS<L1.L)-C2*mC(KL.L)“C1*AS(ICL.L>  72 

41  CONTINUE  73 

C  INTERCHANGE  COLUMNS  74 

DO  44  1=1* N  75 

L=N* 1— I  76 

IF ( InUA (L  *  1 ) -INOX (L * 2) )  42,44,42  77 

42  I  Rw  =  I  Nl)X  ( L  *  1 )  78 

ICL=INUML,2)  79 

DO  4j  k  =  1,n  80 

C1=AC(K.IR*)  81 

C2=AS(K.IRw>  62 

AC(l\*IRrf)=AC(K,ICL>  83 

AS (K, 1K«) =AS <K, ILL)  84 

AC ( K  * ICL ) =C  1  85 

AS ( K • ICL ) =C2  86 

43  CONTINUE  87 

44  CONTINUE  88 

DO  45  i\  =  1,n  89 

IF  ( INOX  (R*  3)  -1 )  4b, 4b * **6  90 

45  CONTINUE  91 

1 0=1  92 

GO  TO  47  93 

46  ID=2  94 

DTC=».0 

DTS=0.0 

47  RETURN  95 

END  96 
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